BOSS 7.1.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtD0ToKSpipipi0pi0.cc
Go to the documentation of this file.
1//------------------------------------------------------------------------
2//
3// Environment:
4// This software is part of the EvtGen package developed jointly
5// for the BaBar and CLEO collaborations. If you use all or part
6// of it, please give an appropriate acknowledgement.
7//
8// Copyright Information: See EvtGen/COPYRIGHT
9// Copyright (C) 1998 Caltech, UCSB
10//
11// Module: EvtD0ToKSpipipi0pi0.cc
12//
13// Description: https://indico.ihep.ac.cn/event/16632/contributions/49342/attachments/23562/26699/meeting.pdf
14//
15// Modification history:
16//
17// Liaoyuan Dong Oct 21, 2022 Module created
18//
19//------------------------------------------------------------------------
20#include "EvtGenBase/EvtPatches.hh"
21#include <fstream>
22#include <stdlib.h>
23#include "EvtGenBase/EvtParticle.hh"
24#include "EvtGenBase/EvtGenKine.hh"
25#include "EvtGenBase/EvtPDL.hh"
26#include "EvtGenBase/EvtReport.hh"
27#include "EvtGenBase/EvtResonance.hh"
28#include "EvtGenBase/EvtResonance2.hh"
29#include "EvtGenModels/EvtD0ToKSpipipi0pi0.hh"
30#include <string>
31#include "EvtGenBase/EvtConst.hh"
32#include "EvtGenBase/EvtDecayTable.hh"
33using std::endl;
34
35EvtD0ToKSpipipi0pi0::~EvtD0ToKSpipipi0pi0() {}
36
37void EvtD0ToKSpipipi0pi0::getName(std::string& model_name){
38 model_name="D0ToKSpipipi0pi0";
39 }
40
41 EvtDecayBase* EvtD0ToKSpipipi0pi0::clone(){
42 return new EvtD0ToKSpipipi0pi0;
43 }
44void EvtD0ToKSpipipi0pi0::init(){
45 checkNArg(0);
46 checkNDaug(5);
47 checkSpinParent(EvtSpinType::SCALAR);
48 checkSpinDaughter(0,EvtSpinType::SCALAR);
49 checkSpinDaughter(1,EvtSpinType::SCALAR);
50 checkSpinDaughter(2,EvtSpinType::SCALAR);
51 checkSpinDaughter(3,EvtSpinType::SCALAR);
52 checkSpinDaughter(4,EvtSpinType::SCALAR);
53
54 string name;
55 Double_t mean_v;
56
57 phi[0] = 0.0;
58 phi[1] = 1.808202793;
59 phi[2] = 1.353154265;
60 phi[3] = 0.02534778524;
61 phi[4] = 2.572695445;
62 phi[5] = 2.572695445;
63 phi[6] = -2.683237531;
64 phi[7] = 0.1209896124;
65
66 rho[0] = 1.0;
67 rho[1] = 0.5632373979;
68 rho[2] = 1.541400902;
69 rho[3] = 163.725261;
70 rho[4] = 3.207856388;
71 rho[5] = 3.207856388;
72 rho[6] = 2.943199635;
73 rho[7] = 3.023783561;
74
75 modetype[0]= 1;
76 modetype[1]= 1;
77 modetype[2]= 1;
78 modetype[3]= 41;
79 modetype[4]= 20;
80 modetype[5]= 33;
81 modetype[6]= 4;
82 modetype[7]= 40;
83
84 //std::cout << "Initializing EvtD0ToKSpipipi0pi0" << std::endl;
85 //for (int i=0; i<8; i++) {
86 // cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
87 //}
88 width1[0] = 0.00849;
89 width1[1] = 0.00849;
90 width1[2] = 0.00849;
91 width1[3] = 0.00849;//ks4pi
92 width1[4] = 0.420;
93 width1[5] = 0.420;
94 width1[6] = 0.142;
95 width1[7] = 0.00849;
96
97 width2[0] = 0.0473;
98 width2[1] = 0.0473;
99 width2[2] = 0.0473;
100 width2[3] = 0.0473;
101 width2[4] = 0.0473;
102 width2[5] = 0.0473;
103 width2[6] = 0.00849;
104 width2[7] = 0.0473;
105
106 mass1[0] = 0.78265;
107 mass1[1] = 0.78265;
108 mass1[2] = 0.78265;
109 mass1[3] = 0.78265;
110 mass1[4] = 1.230;
111 mass1[5] = 1.230;
112 mass1[6] = 1.2295;
113 mass1[7] = 0.78265;
114
115 mass2[0] = 0.89555;
116 mass2[1] = 0.89555;
117 mass2[2] = 0.89555;
118 mass2[3] = 0.89555;
119 mass2[4] = 0.89555;
120 mass2[5] = 0.89555;
121 mass2[6] = 0.78265;
122 mass2[7] = 0.89555;
123
124 mD0M = 1.86486;
125 mD = 1.86486;
126 metap = 0.95778;
127 mkstr = 0.89594;
128 mk0 = 0.497614;
129 mass_Kaon = 0.49368;
130 mass_Pion = 0.13957;
131 mass_Pion2 = 0.0194797849;
132 mass_2Pion = 0.27914;
133 math_2pi = 6.2831852;
134 rD2 = 25.0; // 5*5
135 rDs2 = 25.0; // 5*5
136 rRes2 = 9.0; // 3*3
137 gg1 = 0.5468;
138 gg2 = 0.23; // K*0(1430)
139 GS1 = 0.636619783;
140 GS2 = 0.01860182466;
141 GS3 = 0.1591549458; // 1/(2*math_2pi)
142 GS4 = 0.00620060822;
143
144 int GG[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
145 for (int i=0; i<4; i++) {
146 for (int j=0; j<4; j++) {
147 G[i][j] = GG[i][j];
148 }
149 }
150 double EE[4][4][4][4] =
151{ { {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
152 {{0,0,0,0}, {0,0,0,0}, {0,0,0,1}, {0,0,-1,0}},
153 {{0,0,0,0}, {0,0,0,-1}, {0,0,0,0}, {0,1,0,0} },
154 {{0,0,0,0}, {0,0,1,0}, {0,-1,0,0}, {0,0,0,0} } },
155 { {{0,0,0,0}, {0,0,0,0}, {0,0,0,-1}, {0,0,1,0} },
156 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
157 {{0,0,0,1}, {0,0,0,0}, {0,0,0,0}, {-1,0,0,0}},
158 {{0,0,-1,0}, {0,0,0,0}, {1,0,0,0}, {0,0,0,0} } },
159 { {{0,0,0,0}, {0,0,0,1}, {0,0,0,0}, {0,-1,0,0}},
160 {{0,0,0,-1}, {0,0,0,0}, {0,0,0,0}, {1,0,0,0} },
161 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
162 {{0,1,0,0}, {-1,0,0,0}, {0,0,0,0}, {0,0,0,0} } },
163 { {{0,0,0,0}, {0,0,-1,0}, {0,1,0,0}, {0,0,0,0} },
164 {{0,0,1,0}, {0,0,0,0}, {-1,0,0,0}, {0,0,0,0} },
165 {{0,-1,0,0}, {1,0,0,0}, {0,0,0,0}, {0,0,0,0} },
166 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} } } };
167
168 for (int i=0; i<4; i++) {
169 for (int j=0; j<4; j++) {
170 for (int k=0; k<4; k++) {
171 for (int l=0; l<4; l++) {
172 E[i][j][k][l] = EE[i][j][k][l];
173 }
174 }
175 }
176 }
177
178}
179
180void EvtD0ToKSpipipi0pi0::initProbMax() {
181 setProbMax(52000000);
182}
183
184void EvtD0ToKSpipipi0pi0::decay( EvtParticle *p){
185/*
186 double maxprob = 0.0;
187 for(int ir=0;ir<=60000000;ir++){
188
189 p->initializePhaseSpace(getNDaug(),getDaugs());
190 EvtVector4R D1 = p->getDaug(0)->getP4();
191 EvtVector4R D2 = p->getDaug(1)->getP4();
192 EvtVector4R D3 = p->getDaug(2)->getP4();
193 EvtVector4R D4 = p->getDaug(3)->getP4();
194 EvtVector4R D5 = p->getDaug(4)->getP4();
195
196 double P1[4], P2[4], P3[4], P4[4], P5[4];
197 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
198 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
199 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
200 P4[0] = D4.get(0); P4[1] = D4.get(1); P4[2] = D4.get(2); P4[3] = D4.get(3);
201 P5[0] = D5.get(0); P5[1] = D5.get(1); P5[2] = D5.get(2); P5[3] = D5.get(3);
202
203 double value;
204 int g0[8]={1,1,1,0,1,1,1,1};
205 int g1[8]={1,1,1,0,1,1,1,1};
206 double r0[8]={1,1,1,0,1,1,1,1};
207 double r1[8]={1,1,1,0,1,1,1,1};
208 int spin[8]={0,1,2,0,1020,1020,0,1};
209
210 int nstates=8;
211 calEva(P1, P4, P2, P3, P5, mass1, mass2, width1, width2, rho, phi, r0, r1, g0, g1, spin, modetype, nstates, value);
212
213 if (value<0) continue;
214 if(value>maxprob) {
215 maxprob=value;
216 cout << "ir= " << ir << endl;
217 cout << "double P1[4] = {" << P1[0] <<","<< P1[1] <<","<< P1[2] <<","<< P1[3] <<"};"<< endl;
218 cout << "double P2[4] = {" << P2[0] <<","<< P2[1] <<","<< P2[2] <<","<< P2[3] <<"};"<< endl;
219 cout << "double P3[4] = {" << P3[0] <<","<< P3[1] <<","<< P3[2] <<","<< P3[3] <<"};"<< endl;
220 cout << "double P4[4] = {" << P4[0] <<","<< P4[1] <<","<< P4[2] <<","<< P4[3] <<"};"<< endl;
221 cout << "double P5[4] = {" << P5[0] <<","<< P5[1] <<","<< P5[2] <<","<< P5[3] <<"};"<< endl;
222 cout << "MAX====> " << maxprob << endl;
223
224 }
225 }
226 printf("MAXprob = %.10f\n",maxprob);
227*/
228 p->initializePhaseSpace(getNDaug(),getDaugs());
229 EvtVector4R D1 = p->getDaug(0)->getP4();
230 EvtVector4R D2 = p->getDaug(1)->getP4();
231 EvtVector4R D3 = p->getDaug(2)->getP4();
232 EvtVector4R D4 = p->getDaug(3)->getP4();
233 EvtVector4R D5 = p->getDaug(4)->getP4();
234
235 double P1[4], P2[4], P3[4], P4[4], P5[4];
236 P1[0] = D1.get(0); P1[1] = D1.get(1); P1[2] = D1.get(2); P1[3] = D1.get(3);
237 P2[0] = D2.get(0); P2[1] = D2.get(1); P2[2] = D2.get(2); P2[3] = D2.get(3);
238 P3[0] = D3.get(0); P3[1] = D3.get(1); P3[2] = D3.get(2); P3[3] = D3.get(3);
239 P4[0] = D4.get(0); P4[1] = D4.get(1); P4[2] = D4.get(2); P4[3] = D4.get(3);
240 P5[0] = D5.get(0); P5[1] = D5.get(1); P5[2] = D5.get(2); P5[3] = D5.get(3);
241
242// P1[0] = 0.699073767764678; P1[1] = 0.116528634106469; P1[2] = -0.416631548481574; P1[3] = -0.232030327801847;
243// P2[0] = 0.431675469128871; P2[1] = 0.090021693909732; P2[2] = 0.149419303854623; P2[3] =0.371064840093330;
244// P3[0] = 0.188455625514340; P3[1] = -0.100081780055657;P3[2] = -0.052767368589371; P3[3] = 0.056876884593353;
245// P4[0] = 0.311123417560542; P4[1] = 0.231980078879333; P4[2] = 0.149145989336585; P4[3] = -0.035478344416011;
246// P5[0] = 0.257852771533528; P5[1] = -0.185420043160991;P5[2] = -0.073435020040710; P5[3] = -0.092139879037438;
247
248 double value;
249 int g0[8]={1,1,1,0,1,1,1,1};
250 int g1[8]={1,1,1,0,1,1,1,1};
251 double r0[8]={1,1,1,0,1,1,1,1};
252 double r1[8]={1,1,1,0,1,1,1,1};
253 int spin[8]={0,1,2,0,1020,1020,0,0};
254
255 int nstates=8;
256 calEva(P1, P4, P2, P3, P5, mass1, mass2, width1, width2, rho, phi, r0, r1, g0, g1, spin, modetype, nstates, value);
257 setProb(value);
258
259 return ;
260}
261
262
263void EvtD0ToKSpipipi0pi0::Com_Multi(double a1[2], double a2[2], double res[2])//(a1+b1i)(a2+b2i)
264{
265 res[0] = a1[0]*a2[0]-a1[1]*a2[1];
266 res[1] = a1[1]*a2[0]+a1[0]*a2[1];
267}
268void EvtD0ToKSpipipi0pi0::Com_Divide(double a1[2], double a2[2], double res[2])
269{
270 double tmp = a2[0]*a2[0]+a2[1]*a2[1];
271 res[0] = (a1[0]*a2[0]+a1[1]*a2[1])/tmp;
272 res[1] = (a1[1]*a2[0]-a1[0]*a2[1])/tmp;
273}
274
275//------------base---------------------------------
276double EvtD0ToKSpipipi0pi0::SCADot(double a1[4], double a2[4])
277{
278 double _cal = a1[0]*a2[0]-a1[1]*a2[1]-a1[2]*a2[2]-a1[3]*a2[3];
279 return _cal;
280}
281double EvtD0ToKSpipipi0pi0::Barrier(int l, double sa, double sb, double sc, double r2)
282{
283 double F;
284 double tmp = sa+sb-sc;
285 double q = 0.25*tmp*tmp/sa-sb;
286 if (q < 0) q = -q;
287 double z = q*r2;
288 if (l==1) {
289 F = sqrt(2.0*z/(1.0+z));
290 }
291 else if (l==2) {
292 double z2 = z*z;
293 F = sqrt(13.0*z2/(9.0+3.0*z+z2));
294 } else {
295 F = 1.0;
296 }
297 return F;
298}
299
300double EvtD0ToKSpipipi0pi0::BarrierN(int l, double sa, double sb, double sc, double r, double mass)
301{
302 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
303 if(q < 0) q = -q;
304 double z;
305 //z = q*r*r;
306 z = q*r;
307 double sa0;
308 sa0 = mass*mass;
309 double q0 = (sa0+sb-sc)*(sa0+sb-sc)/(4*sa0)-sb;
310 if(q0 < 0) q0 = -q0;
311 //double z0 = q0*r*r;
312 double z0 = q0*r;
313 double F = 0.0;
314 if(l == 0) F = 1;
315 if(l == 1) F = sqrt((1+z0)/(1+z));
316 if(l == 2) F = sqrt((9+3*z0+z0*z0)/(9+3*z+z*z));
317 return F;
318}
319
320//double EvtD0ToKSpipipi0pi0::BarrierN(int l, double sa, double sb, double sc, double r2)
321//{
322// double B;
323// double tmp = sa + sb -sc;
324// double Q2 = 0.25*tmp*tmp/sa - sb;
325// if (Q2 < 0) Q2 = 1e-16;
326// double Q02 = 0.197321*0.197321/r2;
327// if (l==1){
328// B = sqrt(2.0*Q2/(Q2 + Q02));
329// }
330// else if(l==2){
331// B = sqrt(13.0*Q2*Q2/(Q2*Q2+3.0*Q2*Q02+9.0*Q02*Q02));
332// }
333// else{
334// B = 1.0;
335// }
336// return B;
337//}
338
339//------------------spin-------------------------------------------
340void EvtD0ToKSpipipi0pi0::calt1(double daug1[4], double daug2[4], double t1[4])
341{
342 double p, pq, tmp;
343 double pa[4], qa[4];
344 for(int i=0; i<4; i++) {
345 pa[i] = daug1[i] + daug2[i];
346 qa[i] = daug1[i] - daug2[i];
347 }
348 p = SCADot(pa,pa);
349 pq = SCADot(pa,qa);
350 tmp = pq/p;
351 for(int i=0; i<4; i++) {
352 t1[i] = qa[i] - tmp*pa[i];
353 }
354}
355void EvtD0ToKSpipipi0pi0::calt2(double daug1[4], double daug2[4], double t2[4][4])
356{
357 double p, r;
358 double pa[4], t1[4];
359 calt1(daug1,daug2,t1);
360 r = SCADot(t1,t1)/3.0;
361 for(int i=0; i<4; i++) {
362 pa[i] = daug1[i] + daug2[i];
363 }
364 p = SCADot(pa,pa);
365 for(int i=0; i<4; i++) {
366 for(int j=0; j<4; j++) {
367 t2[i][j] = t1[i]*t1[j] - r*(G[i][j]-pa[i]*pa[j]/p);
368 }
369 }
370}
371void EvtD0ToKSpipipi0pi0::calG2(double daug1[4], double daug2[4], double g2[4][4])
372{
373// double p, r;
374 double p;
375// double pa[4], t1[4];
376 double pa[4];
377// calt1(daug1,daug2,t1);
378// r = SCADot(t1,t1);
379 for(int i=0; i<4; i++)
380 {
381 pa[i] = daug1[i] + daug2[i];
382 }
383 p = SCADot(pa,pa);
384 for(int i=0; i<4; i++)
385 {
386 for(int j=0; j<4; j++)
387 {
388// t2[i][j] = t1[i]*t1[j] - 1.0/3*r*(G[i][j]-pa[i]*pa[j]/p);
389 g2[i][j] = G[i][j]-pa[i]*pa[j]/p;
390 }
391 }
392}
393//-------------------prop--------------------------------------------
394void EvtD0ToKSpipipi0pi0::propagator(double mass2, double mass, double width, double sx, double prop[2])
395{
396 double a[2], b[2];
397 a[0] = 1;
398 a[1] = 0;
399 b[0] = mass2-sx;
400 b[1] = -mass*width;
401 Com_Divide(a,b,prop);
402}
403double EvtD0ToKSpipipi0pi0::wid(double mass2, double mass, double sa, double sb, double sc, double r2, int l)
404{
405 double widm = 0.;
406 double m = sqrt(sa);
407 double tmp = sb-sc;
408 double tmp1 = sa+tmp;
409 double q = 0.25*tmp1*tmp1/sa-sb;
410 if(q<0) q = -q;
411 double tmp2 = mass2+tmp;
412 double q0 = 0.25*tmp2*tmp2/mass2-sb;
413 if(q0<0) q0 = -q0;
414 double z = q*r2;
415 double z0 = q0*r2;
416 double t = q/q0;
417 if(l == 0) {widm = sqrt(t)*mass/m;}
418 else if(l == 1) {widm = t*sqrt(t)*mass/m*(1+z0)/(1+z);}
419 else if(l == 2) {widm = t*t*sqrt(t)*mass/m*(9+3*z0+z0*z0)/(9+3*z+z*z);}
420 return widm;
421}
422double EvtD0ToKSpipipi0pi0::widl1(double mass2, double mass, double sa, double sb, double sc, double r2)
423{
424 double widm = 0.;
425 double m = sqrt(sa);
426 double tmp = sb-sc;
427 double tmp1 = sa+tmp;
428 double q = 0.25*tmp1*tmp1/sa-sb;
429 if(q<0) q = -q;
430 double tmp2 = mass2+tmp;
431 double q0 = 0.25*tmp2*tmp2/mass2-sb;
432 if(q0<0) q0 = -q0;
433 double z = q*r2;
434 double z0 = q0*r2;
435 double F = (1+z0)/(1+z);
436 double t = q/q0;
437 widm = t*sqrt(t)*mass/m*F;
438 return widm;
439}
440void EvtD0ToKSpipipi0pi0::propagatorRBW(double mass2, double mass, double width, double sa, double sb, double sc, double r2, int l, double prop[2])
441{
442 double a[2], b[2];
443 a[0] = 1;
444 a[1] = 0;
445 b[0] = mass2-sa;
446 b[1] = -mass*width*wid(mass2,mass,sa,sb,sc,r2,l);
447 Com_Divide(a,b,prop);
448}
449void EvtD0ToKSpipipi0pi0::propagatorRBWl1(double mass2, double mass, double width, double sa, double sb, double sc, double r2, double prop[2])
450{
451 double a[2], b[2];
452 a[0] = 1;
453 a[1] = 0;
454 b[0] = mass2-sa;
455 b[1] = -mass*width*widl1(mass2,mass,sa,sb,sc,r2);
456 Com_Divide(a,b,prop);
457}
458//------------GS---used by rho----------------------------
459void EvtD0ToKSpipipi0pi0::propagatorGS(double mass2, double mass, double width, double sa, double sb, double sc, double r2, double prop[2])
460{
461 double a[2], b[2];
462 double tmp = sb-sc;
463 double tmp1 = sa+tmp;
464 double q2 = 0.25*tmp1*tmp1/sa-sb;
465 if(q2<0) q2 = -q2;
466
467 double tmp2 = mass2+tmp;
468 double q02 = 0.25*tmp2*tmp2/mass2-sb;
469 if(q02<0) q02 = -q02;
470
471 double q = sqrt(q2);
472 double q0 = sqrt(q02);
473 double m = sqrt(sa);
474 double q03 = q0*q02;
475 double tmp3 = log(mass+2*q0)+1.2760418309; // log(mass_2Pion) = 1.2760418309;
476
477 double h = GS1*q/m*(log(m+2*q)+1.2760418309);
478 double h0 = GS1*q0/mass*tmp3;
479 double dh = h0*(0.125/q02-0.5/mass2)+GS3/mass2;
480 double d = GS2/q02*tmp3+GS3*mass/q0-GS4*mass/q03;
481 double f = mass2/q03*(q2*(h-h0)+(mass2-sa)*q02*dh);
482
483 a[0] = 1.0+d*width/mass;
484 a[1] = 0.0;
485 b[0] = mass2-sa+width*f;
486 b[1] = -mass*width*widl1(mass2,mass,sa,sb,sc,r2);
487 Com_Divide(a,b,prop);
488}
489void EvtD0ToKSpipipi0pi0::propagatorGS_Omg(double mass2, double mass, double width, double sa, double sb, double sc, double r2, double amp_Omg, double phs_Omg, double prop[2])
490{
491 const double Omg_mass = 0.783;
492 const double Omg_width = 0.008;
493 const double Omg_mass2 = 0.613089;
494
495 double tmp_propGS[2];
496 propagatorGS(mass2, mass, width, sa, sb, sc, r2, tmp_propGS);
497 double tmp_propOmg[2];
498 propagatorRBWl1(Omg_mass2, Omg_mass, Omg_width, sa, sb, sc, r2, tmp_propOmg); // Using constant Gamma
499
500 double tmp_propOmg02[2];
501 tmp_propOmg02[0] = Omg_mass2*tmp_propOmg[0];
502 tmp_propOmg02[1] = Omg_mass2*tmp_propOmg[1];
503
504 double tmp_cof[2];
505 tmp_cof[0] = amp_Omg*cos(phs_Omg);
506 tmp_cof[1] = amp_Omg*sin(phs_Omg);
507
508 double amp_M_RBWOm[2];
509 amp_M_RBWOm[0] = tmp_cof[0]*tmp_propOmg02[0]-tmp_cof[1]*tmp_propOmg02[1];
510 amp_M_RBWOm[1] = tmp_cof[1]*tmp_propOmg02[0]+tmp_cof[0]*tmp_propOmg02[1];
511
512 double GS_amp_RBWOm[2];
513 Com_Multi(tmp_propGS, amp_M_RBWOm, GS_amp_RBWOm);
514
515 prop[0] = tmp_propGS[0] + GS_amp_RBWOm[0];
516 prop[1] = tmp_propGS[1] + GS_amp_RBWOm[1];
517
518}
519void EvtD0ToKSpipipi0pi0::rhoab(double sa, double sb, double sc, double res[2]) {
520 double tmp = sa+sb-sc;
521 double q = 0.25*tmp*tmp/sa-sb;
522 if(q>=0) {
523 res[0]=2.0*sqrt(q/sa);
524 res[1]=0.0;
525 } else {
526 res[0]=0.0;
527 res[1]=2.0*sqrt(-q/sa);
528 }
529}
530void EvtD0ToKSpipipi0pi0::rho4Pi(double sa, double res[2]) {
531 double temp = 1.0-0.3116765584/sa; // 0.3116765584=0.13957*0.13957*16
532 if(temp>=0) {
533 res[0]=sqrt(temp)/(1.0+exp(9.8-3.5*sa));
534 res[1]=0.0;
535 } else {
536 res[0]=0.0;
537 res[1]=sqrt(-temp)/(1.0+exp(9.8-3.5*sa));
538 }
539}
540void EvtD0ToKSpipipi0pi0::propagatorsigma500(double sa, double sb, double sc, double prop[2]) {
541 double f = 0.5843+1.6663*sa;
542 const double M = 0.9264;
543 const double mass2 = 0.85821696; // M*M
544 const double mpi2d2 = 0.00973989245;
545 double g1 = f*(sa-mpi2d2)/(mass2-mpi2d2)*exp((mass2-sa)/1.082);
546 double rho1s[2], rho1M[2], rho2s[2], rho2M[2], rho1[2], rho2[2];
547 rhoab(sa,sb,sc,rho1s);
548 rhoab(mass2,sb,sc,rho1M);
549 rho4Pi(sa,rho2s);
550 rho4Pi(mass2,rho2M);
551 Com_Divide(rho1s,rho1M,rho1);
552 Com_Divide(rho2s,rho2M,rho2);
553 double a[2], b[2];
554 a[0] = 1.0;
555 a[1] = 0.0;
556 b[0] = mass2-sa+M*(g1*rho1[1]+0.0024*rho2[1]);
557 b[1] = -M*(g1*rho1[0]+0.0024*rho2[0]);
558 Com_Divide(a,b,prop);
559}
560void EvtD0ToKSpipipi0pi0::KPiSLASS(double sa, double sb, double sc, double prop[2]) {
561/*
562 const double m1430 = 1.463;
563 const double sa0 = 2.140369; // m1430*m1430;
564 const double w1430 = 0.233;
565 const double Lass1 = 0.25/sa0;
566 double tmp = sb-sc;
567 double tmp1 = sa0+tmp;
568 double q0 = Lass1*tmp1*tmp1-sb;
569 if(q0<0) q0 = 1e-16;
570 double tmp2 = sa+tmp;
571 double qs = 0.25*tmp2*tmp2/sa-sb;
572 double q = sqrt(qs);
573 double width = w1430*q*m1430/sqrt(sa*q0);
574 double temp_R = atan(m1430*width/(sa0-sa));
575 if(temp_R<0) temp_R += math_pi;
576 double deltaR = -5.31 + temp_R;
577 double temp_F = atan(2.14*q/(2.0-1.926*qs)); // 2.0*1.07 = 2.14; 1.8*1.07 = 1.926
578 if(temp_F<0) temp_F += math_pi;
579 double deltaF = 2.33 + temp_F;
580 double deltaS = deltaR + 2.0*deltaF;
581 double t1 = 0.8*sin(deltaF);
582 double t2 = sin(deltaR);
583 double CF[2], CS[2];
584 CF[0] = cos(deltaF);
585 CF[1] = sin(deltaF);
586 CS[0] = cos(deltaS);
587 CS[1] = sin(deltaS);
588 prop[0] = t1*CF[0] + t2*CS[0];
589 prop[1] = t1*CF[1] + t2*CS[1];
590*/
591 const double m1430 = 1.441;
592 const double sa0 = 1.441*1.441; // m1430*m1430;
593 const double w1430 = 0.193;
594 const double Lass1 = 0.25/sa0;
595 double tmp = sb-sc;
596 double tmp1 = sa0+tmp;
597 double q0 = Lass1*tmp1*tmp1-sb;
598 if(q0<0) q0 = -q0;
599 double tmp2 = sa+tmp;
600 double qs = 0.25*tmp2*tmp2/sa-sb;
601 double q = sqrt(qs);
602 double width = w1430*q*m1430/sqrt(sa*q0);
603 double temp_R = atan(m1430*width/(sa0-sa));
604 if(temp_R<0) temp_R += math_pi;
605 double deltaR = -1.915 + temp_R; //fiR=-109.7
606 double temp_F = atan(0.226*q/(2.0-3.819*qs)); // 2.0*0.113 = 0.226; 0.113*33.8 = 3.819
607 if(temp_F<0) temp_F += math_pi;
608 double deltaF = 0.002 + temp_F; //fiF=0.1
609 double deltaS = deltaR + 2.0*deltaF;
610 double t1 = 0.96*sin(deltaF);
611 double t2 = sin(deltaR);
612 double CF[2], CS[2];
613 CF[0] = cos(deltaF);
614 CF[1] = sin(deltaF);
615 CS[0] = cos(deltaS);
616 CS[1] = sin(deltaS);
617 prop[0] = t1*CF[0] + t2*CS[0];
618 prop[1] = t1*CF[1] + t2*CS[1];
619
620}
621
622//0-15: 4P*4;
623//16-23: non-resonance;
624//24: pdf after calculating;
625//28: mass and width;
626//29-30: amp, phase;
627//31-34: g0: do/not cal prop about KaPi1; g1: do/not cal prop about Pi2Pim; modetype:
628//35-38: r0/1: centrifugal barrier; nstates: number of state;
629//39-40: use [first:last] amplitudes.
630
631void EvtD0ToKSpipipi0pi0::calEva( double* PKs,double* Pip,double* Pim,double* Pi01,double* Pi02,double* mass1, double* mass2, double* width1, double* width2, double* amp, double* phase, double* r0, double* r1, int* g0, int* g1, int* g2, int* modetype,double nstates,double &Result)
632{
633
634
635 //Generally, M(K*0(1430)) is been set to a constant(result from CLEO)
636 //modetype: 0, a1(1260)K.
637 // 1, K1(1270)Pi and K1(1270)->K*Pi
638 // 2, VV(2,10,11,12)
639 // 3, K1(1270)Pi and K1(1270)->rhoK(3,14)
640 // 4, K1(1270)Pi and K1(1270)->K*0(1430)Pi
641 // 20, (KPi)Srho(20,22)
642 // 21, K*(PiPi)S(21,32)
643 // 42, (KPi)S(PiPi)S
644 // 13, (K*Pi)_{P}pi
645 // 5, (rhoPi)_{P}K
646 // 6, (rhoPi)_{V}K
647 // 7, (K*Pi)_{V}Pi
648
649 double cof[2], amp_tmp[2],amp_tmp1[2], amp_tmp2[2], amp_PDF[2], PDF[2];
650 double pomega1[4], pomega2[4], pkst1[4], pkst2[4],pD0[4],prho0[4],prhom1[4],prhom2[4],prhop1[4],prhop2[4],pb10[4],pa101[4],pa102[4],Pf01[4],Pf02[4],pk14001[4],pk14002[4],pk1270[4],pk1400[4],pf1285[4],prho01[4],prho02[4],pkstm[4],pkpi01SW[4],pkpi02SW[4],p1omega1[4],p1omega2[4],prho1450[4],pa1p[4],pkstm1[4],pkstm2[4],pk127m1[4],pk127m2[4],pk127p1[4],pk127p2[4],prhom11[4],prhom12[4],prhom111[4],prhom112[4];
651
652 for(int i=0;i!=4;i++) {
653
654 pomega1[i] = Pi01[i]+Pip[i] + Pim[i];
655 pkst2[i] = PKs[i] + Pi02[i];
656 pomega2[i] =Pip[i] + Pim[i] + Pi02[i];
657 pkst1[i] = PKs[i] + Pi01[i];
658 pD0[i] = PKs[i] + Pi01[i] + Pip[i] + Pim[i] + Pi02[i];
659 prho0[i] = Pip[i] + Pim[i];
660 prho01[i] = Pip[i] + Pim[i];
661 prho02[i] = Pip[i] + Pim[i];
662 prhom1[i] = Pi01[i] + Pim[i];
663 prhom2[i] = Pim[i] + Pi02[i];
664 prhop1[i] = Pi01[i] + Pip[i];
665 prhop2[i] = Pip[i] + Pi02[i];
666 pb10[i] = Pi01[i] + Pip[i] + Pim[i] + Pi02[i];
667
668 //a1p kstm
669 pkstm1[i] = PKs[i] + Pim[i];
670 pkstm2[i] = PKs[i] + Pim[i];
671 pkstm[i] = PKs[i] + Pim[i];
672 //a1m kstp
673 //a10 kst0
674 pa101[i] =Pi01[i] + Pip[i] + Pim[i] ;
675 pa102[i] = Pip[i] + Pim[i] + Pi02[i];
676 pa1p[i] = Pi01[i] + Pip[i] + Pi02[i];
677 //f0500
678 Pf01[i]=Pi01[i]+Pi02[i];
679 Pf02[i]=Pip[i]+Pim[i];
680 //k1(1270)
681 //k1(1400)
682 pk14001[i] = PKs[i] + Pi01[i] + Pip[i] + Pim[i];
683 pk14002[i] = PKs[i] + Pi02[i] + Pip[i] + Pim[i];
684 //k1(1270)->KS rho
685 pk1270[i]=PKs[i] + Pip[i] + Pim[i];
686 //k1(1400)->kstar pi0
687 pk1400[i] = PKs[i] + Pi01[i] + Pi02[i];
688 //f1(1285)->4pi
689 pf1285[i] = Pi01[i] + Pip[i] + Pim[i] + Pi02[i];
690 // omega KPi01SW
691 p1omega1[i] = Pip[i] + Pim[i] + Pi01[i];
692 p1omega2[i] = Pip[i] + Pim[i] + Pi02[i];
693 pkpi01SW[i] = PKs[i] + Pi01[i];
694 pkpi02SW[i] = PKs[i] + Pi02[i];
695 // prho1450
696 prho1450[i] = Pi01[i] + Pip[i] + Pim[i] + Pi02[i];
697 //pk12701 pk12702
698 pk127m1[i] = PKs[i] + Pi01[i] + Pim[i];
699 pk127m2[i] = PKs[i] + Pi02[i] + Pim[i];
700 pk127p1[i] = PKs[i] + Pi01[i] + Pip[i];
701 pk127p2[i] = PKs[i] + Pi02[i] + Pip[i];
702 prhom11[i] = Pi01[i] + Pim[i];
703 prhom12[i] =Pim[i] + Pi02[i];
704 prhom111[i] = Pi01[i] + Pim[i];
705 prhom112[i] =Pim[i] + Pi02[i];
706
707 }
708
709 double sKs,sPi01,skstm,sPim,sPip,sPi02,somega1,sk1270,somega2,skst1,skst2,sD0,srho0,srhop1,srhop2,srhom1,srhom2,sb10,skstm1,skstm2,sa1p,sa101,sa102,sf01,sf02,sk14001,sk14002,sk1400,sf1285,srho01,srho02,s1omega1,s1omega2,skpi01SW,skpi02SW,srho1450,sk127p1,sk127p2,sk127m1,sk127m2,srhom11,srhom12,srhom111,srhom112;
710
711 sKs = SCADot(PKs,PKs);
712 sPi01 = SCADot(Pi01,Pi01);
713 sPip = SCADot(Pip,Pip);
714 sPim = SCADot(Pim,Pim);
715 sPi02 = SCADot(Pi02,Pi02);
716 somega1 = SCADot(pomega1,pomega1);
717 somega2 = SCADot(pomega2,pomega2);
718 skst1 = SCADot(pkst1,pkst1);
719 skst2 = SCADot(pkst2,pkst2);
720 sD0 = SCADot(pD0,pD0);
721 srhop1 =SCADot(prhop1,prhop1);
722 srhop2 =SCADot(prhop2,prhop2);
723 srhom1 =SCADot(prhom1,prhom1);
724 srhom2 =SCADot(prhom2,prhom2);
725 srho0 =SCADot(prho0,prho0);
726 srho01 =SCADot(prho01,prho01);
727 srho02 =SCADot(prho02,prho02);
728 sb10 =SCADot(pb10,pb10);
729 //a1p kstm
730 skstm1 = SCADot(pkstm1,pkstm1);
731 skstm2 = SCADot(pkstm2,pkstm2);
732 skstm = SCADot(pkstm,pkstm);
733 sk1270 = SCADot(pk1270,pk1270);
734 sa1p = SCADot(pa1p,pa1p);
735 //a1m kstp
736 //a10 kst0
737 sa101 = SCADot(pa101,pa101);
738 sa102 = SCADot(pa102,pa102);
739 // sf01 sf02
740 sf01 = SCADot(Pf01,Pf01);
741 sf02 = SCADot(Pf02,Pf02);
742 //sk1p sk2p
743 //sk14001 sk14002
744 sk14001 =SCADot(pk14001,pk14001);
745 sk14002 =SCADot(pk14002,pk14002);
746
747 //sk1400
748 sk1400=SCADot(pk1400,pk1400);
749 //sf1285
750 sf1285=SCADot(pf1285,pf1285);
751 // omegga KPi01SW
752 s1omega1=SCADot(p1omega1,p1omega1);
753 s1omega2=SCADot(p1omega2,p1omega2);
754 skpi01SW=SCADot(pkpi01SW,pkpi01SW);
755 skpi02SW=SCADot(pkpi02SW,pkpi02SW);
756 //srho1450
757 srho1450=SCADot(prho1450,prho1450);
758 ///sk1270 to k*0 pi-
759 sk127p1 = SCADot(pk127p1,pk127p1);
760 sk127p2 = SCADot(pk127p2,pk127p2);
761 sk127m1 = SCADot(pk127m1,pk127m1);
762 sk127m2 = SCADot(pk127m2,pk127m2);
763 srhom11 = SCADot(prhom11,prhom11);
764 srhom12 = SCADot(prhom12,prhom12);
765 srhom111 = SCADot(prhom111,prhom111);
766 srhom112 = SCADot(prhom112,prhom112);
767
768
769////////////////////////////////////////////////////////////////
770 double t1rho01[4],t1rho02[4],t1rhop1[4],t1rhop2[4],t1rhom1[4],t1rhom2[4],t1omega01[4],t1omega02[4],t1omegap1[4],t1omegap2[4],t1omegam1[4],t1omegam2[4],t1kst1[4],t1kst2[4],t1D0_b1[4],t1D0_k14001[4],t1D0_k14002[4],t1D0[4],t1D0_k1400[4],t1D0KPi01SW[4],t1D0KPi02SW[4],t1rho0[4],t1D0_rho1450ks0[4],t1rho14501[4],t1rho14502[4],t1D0_f1285[4],t1f1285[4],t1rhom11[4],t1rhom12[4],t1rhom111[4],t1rhom112[4];
771 double p1omega01[4][4],p1omega02[4][4],p1omegap1[4][4],p1omegap2[4][4],p1omegam1[4][4],p1omegam2[4][4];
772 double p1b1_1[4][4],p1b1_2[4][4],pk1400_1[4][4],pk1400_2[4][4],p2a1270[4][4],p1rho14501[4][4],p1rho14502[4][4],p1a1m_1[4][4],p1a1m_2[4][4],t2_a102[4][4];
773 double p2k1400_1[4][4],p2k1400_2[4][4],p2f1285[4][4],t2f1285[4][4];
774 //Sn---L(i)
775 calt1(Pip,Pim,t1rho01); calt1(Pip,Pim,t1rho02);calt1(Pip,Pim,t1rho0);
776 calt1(Pip,Pi01,t1rhop1); calt1(Pip,Pi02,t1rhop2);
777 calt1(Pi01,Pim,t1rhom1);calt1(Pi02,Pim,t1rhom2);
778 calt1(prho0,Pi01,t1omega01);calt1(prho0,Pi02,t1omega02);
779 calt1(prhop1,Pim,t1omegap1);calt1(prhop2,Pim,t1omegap2);
780 calt1(prhom1,Pip,t1omegam1);calt1(prhom2,Pip,t1omegam2);
781 calt1(PKs,Pi01,t1kst1); calt1(PKs,Pi02,t1kst2);
782 calt1(p1omega1,pkpi02SW,t1D0KPi02SW);
783 calt1(p1omega2,pkpi01SW,t1D0KPi01SW);
784 calt1(pb10,PKs,t1D0_b1);
785 calt1(pk14001,Pi02,t1D0_k14001);
786 calt1(pk14002,Pi01,t1D0_k14002);
787
788 //Sn---P(ij)
789 calG2(prho0,Pi01,p1omega01);calG2(prho0,Pi02,p1omega02);calG2(prhop1,Pim,p1omegap1);
790 calG2(prhop2,Pim,p1omegap2);calG2(prhom1,Pip,p1omegam1);calG2(prhom2,Pip,p1omegam2);
791
792 calG2(pomega1,Pi02,p1b1_1); calG2(pomega2,Pi01,p1b1_2);
793 calG2(pomega1,PKs,pk1400_1); calG2(pomega2,PKs,pk1400_2);
794 // calG2(pkst1,Pi02,p2k1400_1); calG2(pkst2,Pi01,p2k1400_2);
795
796 PDF[0] = 0.0; PDF[1] = 0.0;
797
798 ////
799 double t1D0_omega1kst2[4],t1D0_omega2kst1[4];
800 calt1(pomega1,pkst2,t1D0_omega1kst2); calt1(pomega2,pkst1,t1D0_omega2kst1);
801
802 //Sn(l=2)
803 double t2D0_omega1kst2[4][4],t2D0_omega2kst1[4][4],t2b1_1[4][4],t2b1_2[4][4],t2k1400_1[4][4],t2k1400_2[4][4],t2k1270[4][4],t2D0_k1400[4][4];
804 calt2(pomega1,pkst2,t2D0_omega1kst2); calt2(pomega2,pkst1,t2D0_omega2kst1);
805 calt2(pomega1,Pi02,t2b1_1); calt2(pomega2,Pi01,t2b1_2);
806 calt2(pomega1,PKs,t2k1400_1); calt2(pomega2,PKs,t2k1400_2);
807
808 // rho from omega
809 double pro1V0[2],pro1Vm2[2],pro1Vp1[2],pro1Vp2[2],pro1Vm1[2],pro1V01[2],pro1V02[2],pro11Vp1[2],pro11Vp2[2],pro11Vm1[2],pro11Vm2[2],pro111Vm1[2],pro111Vm2[2];
810
811 propagatorGS(0.60102807,0.77526,0.1478,srho0,sPip,sPim,rRes2,pro1V0);
812 propagatorGS(0.60102807,0.77526,0.1478,srho01,sPip,sPim,rRes2,pro1V01);
813 propagatorGS(0.60102807,0.77526,0.1478,srho02,sPip,sPim,rRes2,pro1V02);
814 propagatorGS(0.60102807,0.77526,0.1478,srhop1,sPip,sPi01,rRes2,pro1Vp1);
815 propagatorGS(0.60102807,0.77526,0.1478,srhop2,sPip,sPi02,rRes2,pro1Vp2);
816 propagatorGS(0.60102807,0.77526,0.1478,srhom1,sPi01,sPim,rRes2,pro1Vm1);
817 propagatorGS(0.60102807,0.77526,0.1478,srhom2,sPi02,sPim,rRes2,pro1Vm2);
818 propagatorGS(0.60102807,0.77526,0.1478,srhom1,sPim,sPi01,rRes2,pro11Vm1);
819 propagatorGS(0.60102807,0.77526,0.1478,srhom2,sPim,sPi02,rRes2,pro11Vm2);
820
821
822 double B1_V01(-1),B1_V02(-1),B1_Vp1(-1),B1_Vp2(-1),B1_Vm1(-1),B1_Vm2(-1),B1_V0(-1),B1_Vm111(-1),B1_Vm112(-1),B1_Vm11(-1),B1_Vm12(-1);
823 // rho from omega Blatt-We
824 B1_V01=BarrierN(1,srho01,sPip,sPim,rRes2,0.77526);
825 B1_V02=BarrierN(1,srho02,sPip,sPim,rRes2,0.77526);
826 B1_V0=BarrierN(1,srho0,sPip,sPim,rRes2,0.77526);
827 B1_Vp1=BarrierN(1,srhop1,sPip,sPi01,rRes2,0.77526);
828 B1_Vp2=BarrierN(1,srhop2,sPip,sPi02,rRes2,0.77526);
829 B1_Vm1=BarrierN(1,srhom1,sPi01,sPim,rRes2,0.77526);
830 B1_Vm2=BarrierN(1,srhom2,sPi02,sPim,rRes2,0.77526);
831 B1_Vm111=BarrierN(1,srhom111,sPim,sPi01,rRes2,0.77526);
832 B1_Vm112=BarrierN(1,srhom112,sPim,sPi02,rRes2,0.77526);
833 B1_Vm11=BarrierN(1,srhom11,sPim,sPi01,rRes2,0.77526);
834 B1_Vm12=BarrierN(1,srhom12,sPim,sPi02,rRes2,0.77526);
835
836
837 double B1_omega10(-1),B1_omega1p(-1),B1_omega1m(-1),B1_omega20(-1),B1_omega2p(-1),B1_omega2m(-1),B1_D0_k1400(-1),B1_D0_k14001(-1),B1_D0_k14002(-1),B2_k1400_1(-1),B2_k1400(-1),B2_k1400_2(-1),B1_D0_k1270(-1),B2_k1270(-1),B1_D0_k1400rho0(-1),B2_D0_k1400rho0(-1),B1_D0_f1(-1),B1_D0_a1pkpim1(-1),B1_D0_a1pkpim2(-1),B1_f1285(-1);
838
839 double B1_D0_omega1(-1),B1_D0_omega2(-1),B1_kst1(-1),B1_kst2(-1),B1_D0_b1(-1),B1_a11p(-1),B1_a11m(-1),B2_a12p(-1),B2_a12m(-1),B1_kstm1(-1),B1_kstm2(-1);
840
841
842 double B2_D0_omega1(-1),B2_D0_omega2(-1),B2_b1_1(-1),B2_b1_2(-1);
843
844 double temp_PDF11, temp_PDF12, temp_PDF13, temp_PDF21, temp_PDF22, temp_PDF23;
845 double pro0_01[2],pro0_p1[2],pro0_m1[2],pro0_02[2],pro0_p2[2],pro0_m2[2];
846
847 double proa_p1[2],proa_m1[2],proa_p2[2],proa_m2[2];
848 double pro1_1[2],pro1_2[2],pro1_3[2],pro1_4[2],pf0_1[2],pf0_2[2],pro00_01[2],pro00_02[2];
849
850 //double pro0_1[2],pro0_2[2],pro0_3[2],pro0_4[2];
851 double pro1_01[2],pro1_p1[2],pro1_m1[2],pro1_02[2],pro1_p2[2],pro1_m2[2];
852 double prob1_1[2],prob1_2[2];
853
854 double tpro_01[2],tpro_p1[2],tpro_m1[2],tpro_02[2],tpro_p2[2],tpro_m2[2];
855
856 double pro_11[2],pro_12[2],pro_13[2],pro_21[2],pro_22[2],pro_23[2],pro_a1p[2],pro_a101[2],pro_a102[2],pfo_a1p[2],pfo_a101[2],pfo_a102[2];
857 double mass1sq, mass2sq;
858
859
860
861 double B2_a1p_rhop1(-1),B1_D0_a1pkstm(-1),B1_D0_a1pkstm1(-1),B1_D0_a1pkstm2(-1),B2_a1p_rhop2(-1),B2_D0_a1pkstm(-1),B1_kst(-1),B1_kstm(-1),B1_kst01(-1),B1_kst02(-1),B1_rho14501(-1),B1_rho14502(-1);
862 double B1_a10_f0(-1),B1_a101(-1),B1_a102(-1),B1_D0_a1kst(-1);
863 double B1_a10_kst0(-1),B1_a10_kst(-1),B2_a10_rho(-1),B1_D0_a10kst0(-1),B2_D0_a10kst0(-1),B1_D0_a10kst01(-1),B1_D0_a10kst02(-1),B2_D0_a10kst01(-1),B2_D0_a10kst02(-1);
864 double p1a1p_1[4][4],p1a1p_2[4][4],t2_a1p1[4][4],t2_a1p2[4][4],t2D0_a1p[4][4],p1_a102[4][4],p1_a103[4][4],p2_a102[4][4],p2_a103[4][4],p1_k12701[4][4],p1_k12702[4][4],p1_k127p1[4][4],p1_k127p2[4][4],p1_k127m1[4][4],p1_k127m2[4][4];
865 double t2D0_a101[4][4],t2D0_a102[4][4];
866 double t1_a102[4][4],t1_a103[4][4],t2_a103[4][4];
867 double t1kstm1[4],t1kstm2[4],t1kstm[4],t1D0_a1p[4],t1D0_a1p1[4],t1D0_a1p2[4],t2D0a101[4][4];
868 double t1a1p[4],t1a101[1],t1a102[1],t1D0_kpimSW1[4],t1D0_kpimSW2[4];
869 double t1D0_a101[4],t1D0_a102[4],t2D0a102[4][4],t1D0a101[4],t1D0a102[4];
870 double temp_PDF1,temp_PDF2;
871
872
873 //---------------------------------------------------------------------------
874
875 for(int i=0; i<nstates; i++) {
876 cof[0] = amp[i]*cos(phase[i]);
877 cof[1] = amp[i]*sin(phase[i]);
878 mass1sq = mass1[i]*mass1[i];
879 mass2sq = mass2[i]*mass2[i];
880 temp_PDF11 = 0;
881 temp_PDF12 = 0;
882 temp_PDF13 = 0;
883 temp_PDF21 = 0;
884 temp_PDF22 = 0;
885 temp_PDF23 = 0;
886 temp_PDF1 = 0;
887 temp_PDF2 = 0;
888
889 //D to V1V2, V1--PKsPi01, V2--V3Pip, V4--PimPi02
890 if(modetype[i]==1){
891 if (g2[i]==0)//spin=0 D to rho omega
892 {
893 if(B1_omega10==-1) {
894
895 B1_omega10=BarrierN(1,somega1,srho01,sPi01,rRes2,mass1[i]);
896 B1_omega1p=BarrierN(1,somega1,srhop1,sPim,rRes2,mass1[i]);
897 B1_omega1m=BarrierN(1,somega1,srhom1,sPip,rRes2,mass1[i]);
898
899 B1_omega20=BarrierN(1,somega2,srho02,sPi02,rRes2,mass1[i]);
900 B1_omega2p=BarrierN(1,somega2,srhop2,sPim,rRes2,mass1[i]);
901 B1_omega2m=BarrierN(1,somega2,srhom2,sPip,rRes2,mass1[i]);
902 }
903 if(B1_kst1==-1) {
904 B1_kst1=BarrierN(1,skst1,sKs,sPi01,rRes2,mass2[i]);
905 B1_kst2=BarrierN(1,skst2,sKs,sPi02,rRes2,mass2[i]);
906 }
907 for(int i=0; i<4; i++)
908 {
909 for(int j=0; j<4; j++)
910 {
911 double tempomega1 = pomega1[j]*G[j][j]*G[i][i];// p(omega)
912 double tempomega2 = pomega2[j]*G[j][j]*G[i][i];
913
914 for(int k=0; k<4; k++)
915 {
916 double temt1omega01 = t1omega01[k]*G[k][k];//Sn(ome)
917 double temt1omega02 = t1omega02[k]*G[k][k];
918 double temt1omegap1 = t1omegap1[k]*G[k][k];
919 double temt1omegap2 = t1omegap2[k]*G[k][k];
920 double temt1omegam1 = t1omegam1[k]*G[k][k];
921 double temt1omegam2 = t1omegam2[k]*G[k][k];
922
923 for(int l=0; l<4; l++)
924 {
925 if((i==j)||(i==k)||(i==l)||(j==k)||(j==l)||(k==l)) continue;
926
927 double temt1rho01 = t1rho01[l]*G[l][l]*E[i][j][k][l];
928 double temt1rho02 = t1rho02[l]*G[l][l]*E[i][j][k][l];
929 double temt1rhop1 = t1rhop1[l]*G[l][l]*E[i][j][k][l];
930 double temt1rhop2 = t1rhop2[l]*G[l][l]*E[i][j][k][l];
931 double temt1rhom1 = t1rhom1[l]*G[l][l]*E[i][j][k][l];
932 double temt1rhom2 = t1rhom2[l]*G[l][l]*E[i][j][k][l];
933
934 for(int m=0; m<4; m++)
935 {
936 temp_PDF11 += G[m][m]*t1kst1[m]*p1omega02[m][i]*tempomega2*temt1omega02*temt1rho01;
937 temp_PDF12 += G[m][m]*t1kst1[m]*p1omegap2[m][i]*tempomega2*temt1omegap2*temt1rhop2;
938 temp_PDF13 += G[m][m]*t1kst1[m]*p1omegam2[m][i]*tempomega2*temt1omegam2*temt1rhom2;
939
940 temp_PDF21 += G[m][m]*t1kst2[m]*p1omega01[m][i]*tempomega1*temt1omega01*temt1rho02;
941 temp_PDF22 += G[m][m]*t1kst2[m]*p1omegap1[m][i]*tempomega1*temt1omegap1*temt1rhop1;
942 temp_PDF23 += G[m][m]*t1kst2[m]*p1omegam1[m][i]*tempomega1*temt1omegam1*temt1rhom1;
943 }
944 }
945 }
946 }
947 }
948
949 temp_PDF11 = temp_PDF11*B1_omega20*B1_kst1*B1_V02;
950 temp_PDF12 = temp_PDF12*B1_omega2p*B1_kst1*B1_Vp2;
951 temp_PDF13 = temp_PDF13*B1_omega2m*B1_kst1*B1_Vm2;
952 temp_PDF21 = temp_PDF21*B1_omega10*B1_kst2*B1_V01;
953 temp_PDF22 = temp_PDF22*B1_omega1p*B1_kst2*B1_Vp1;
954 temp_PDF23 = temp_PDF23*B1_omega1m*B1_kst2*B1_Vm1;
955 }
956 else if (g2[i]==1)
957 {
958
959 if(B1_omega10==-1) {
960
961 B1_omega10=BarrierN(1,somega1,srho01,sPi01,rRes2,mass1[i]);
962 B1_omega1p=BarrierN(1,somega1,srhop1,sPim,rRes2,mass1[i]);
963 B1_omega1m=BarrierN(1,somega1,srhom1,sPip,rRes2,mass1[i]);
964
965 B1_omega20=BarrierN(1,somega2,srho02,sPi02,rRes2,mass1[i]);
966 B1_omega2p=BarrierN(1,somega2,srhop2,sPim,rRes2,mass1[i]);
967 B1_omega2m=BarrierN(1,somega2,srhom2,sPip,rRes2,mass1[i]);
968 }
969
970 if(B1_kst1==-1) {
971
972 B1_kst1=BarrierN(1,skst1,sKs,sPi01,rRes2,mass2[i]);
973 B1_kst2=BarrierN(1,skst2,sKs,sPi02,rRes2,mass2[i]);
974 }
975
976 if(B1_D0_omega1==-1) {
977
978 B1_D0_omega1=BarrierN(1,sD0,somega1,skst2,rDs2,mD0M);
979 B1_D0_omega2=BarrierN(1,sD0,somega2,skst1,rDs2,mD0M);
980 }
981 for(int i=0; i<4; i++)
982 {
983 double tempD0 = pD0[i]*G[i][i];
984
985 for(int j=0; j<4; j++)
986 {
987 double temt1kst1 = t1kst1[j]*G[j][j];
988 double temt1kst2 = t1kst2[j]*G[j][j];
989
990 for(int k=0; k<4; k++)
991 {
992 double temt1D0_omega1kst2 = t1D0_omega1kst2[k]*G[k][k];
993 double temt1D0_omega2kst1 = t1D0_omega2kst1[k]*G[k][k];
994
995 for(int l=0; l<4; l++)
996 {
997 for(int m=0; m<4; m++)
998 {
999 double temp1omega01 = p1omega01[l][m]*G[l][l]*G[m][m]*E[i][j][k][l];
1000 double temp1omega02 = p1omega02[l][m]*G[l][l]*G[m][m]*E[i][j][k][l];
1001 double temp1omegap1 = p1omegap1[l][m]*G[l][l]*G[m][m]*E[i][j][k][l];
1002 double temp1omegap2 = p1omegap2[l][m]*G[l][l]*G[m][m]*E[i][j][k][l];
1003 double temp1omegam1 = p1omegam1[l][m]*G[l][l]*G[m][m]*E[i][j][k][l];
1004 double temp1omegam2 = p1omegam2[l][m]*G[l][l]*G[m][m]*E[i][j][k][l];
1005
1006 for(int n=0; n<4; n++)
1007 {
1008 double tempomega1 = pomega1[n]*G[n][n];
1009 double tempomega2 = pomega2[n]*G[n][n];
1010
1011 for(int p=0; p<4; p++)
1012 {
1013 double temt1omega01 = t1omega01[p]*G[p][p];
1014 double temt1omega02 = t1omega02[p]*G[p][p];
1015 double temt1omegap1 = t1omegap1[p]*G[p][p];
1016 double temt1omegap2 = t1omegap2[p]*G[p][p];
1017 double temt1omegam1 = t1omegam1[p]*G[p][p];
1018 double temt1omegam2 = t1omegam2[p]*G[p][p];
1019
1020 for(int q=0; q<4; q++)
1021 {
1022 if((i==j)||(i==k)||(i==l)||(j==k)||(j==l)||(k==l)) continue;
1023 if((m==n)||(m==p)||(m==q)||(n==p)||(n==q)||(p==q)) continue;
1024
1025 temp_PDF11 += tempD0*temt1kst1*temt1D0_omega2kst1*temp1omega02*E[m][n][p][q]*tempomega2*temt1omega02*t1rho01[q]*G[q][q];
1026 temp_PDF12 += tempD0*temt1kst1*temt1D0_omega2kst1*temp1omegap2*E[m][n][p][q]*tempomega2*temt1omegap2*t1rhop2[q]*G[q][q];
1027 temp_PDF13 += tempD0*temt1kst1*temt1D0_omega2kst1*temp1omegam2*E[m][n][p][q]*tempomega2*temt1omegam2*t1rhom2[q]*G[q][q];
1028 temp_PDF21 += tempD0*temt1kst2*temt1D0_omega1kst2*temp1omega01*E[m][n][p][q]*tempomega1*temt1omega01*t1rho02[q]*G[q][q];
1029 temp_PDF22 += tempD0*temt1kst2*temt1D0_omega1kst2*temp1omegap1*E[m][n][p][q]*tempomega1*temt1omegap1*t1rhop1[q]*G[q][q];
1030 temp_PDF23 += tempD0*temt1kst2*temt1D0_omega1kst2*temp1omegam1*E[m][n][p][q]*tempomega1*temt1omegam1*t1rhom1[q]*G[q][q];
1031
1032 }
1033 }
1034 }
1035 }
1036 }
1037 }
1038 }
1039 }
1040 temp_PDF11 = temp_PDF11*B1_D0_omega2*B1_omega20*B1_kst1*B1_V02;
1041 temp_PDF12 = temp_PDF12*B1_D0_omega2*B1_omega2p*B1_kst1*B1_Vp2;
1042 temp_PDF13 = temp_PDF13*B1_D0_omega2*B1_omega2m*B1_kst1*B1_Vm2;
1043 temp_PDF21 = temp_PDF21*B1_D0_omega1*B1_omega10*B1_kst2*B1_V01;
1044 temp_PDF22 = temp_PDF22*B1_D0_omega1*B1_omega1p*B1_kst2*B1_Vp1;
1045 temp_PDF23 = temp_PDF23*B1_D0_omega1*B1_omega1m*B1_kst2*B1_Vm1;
1046 }
1047 else if (g2[i]==2)
1048 {
1049 if(B1_omega10==-1) {
1050 B1_omega10=BarrierN(1,somega1,srho01,sPi01,rRes2,mass1[i]);
1051 B1_omega1p=BarrierN(1,somega1,srhop1,sPim,rRes2,mass1[i]);
1052 B1_omega1m=BarrierN(1,somega1,srhom1,sPip,rRes2,mass1[i]);
1053
1054 B1_omega20=BarrierN(1,somega2,srho02,sPi02,rRes2,mass1[i]);
1055 B1_omega2p=BarrierN(1,somega2,srhop2,sPim,rRes2,mass1[i]);
1056 B1_omega2m=BarrierN(1,somega2,srhom2,sPip,rRes2,mass1[i]);
1057
1058 }
1059
1060 if(B1_kst1==-1) {
1061 B1_kst1=BarrierN(1,skst1,sKs,sPi01,rRes2,mass2[i]);
1062 B1_kst2=BarrierN(1,skst2,sKs,sPi02,rRes2,mass2[i]);
1063 }
1064
1065 if(B2_D0_omega1==-1){
1066 B2_D0_omega1=BarrierN(2,sD0,somega1,skst2,rDs2,mD0M);
1067 B2_D0_omega2=BarrierN(2,sD0,somega2,skst1,rDs2,mD0M);
1068 }
1069 for(int i=0; i<4; i++)
1070 {
1071 for(int j=0; j<4; j++)
1072 {
1073 double tempomega1 = pomega1[j]*G[j][j]*G[i][i];
1074 double tempomega2 = pomega2[j]*G[j][j]*G[i][i];
1075
1076 for(int k=0; k<4; k++)
1077 {
1078 double temt1omega01 = t1omega01[k]*G[k][k];
1079 double temt1omega02 = t1omega02[k]*G[k][k];
1080 double temt1omegap1 = t1omegap1[k]*G[k][k];
1081 double temt1omegap2 = t1omegap2[k]*G[k][k];
1082 double temt1omegam1 = t1omegam1[k]*G[k][k];
1083 double temt1omegam2 = t1omegam2[k]*G[k][k];
1084
1085 for(int l=0; l<4; l++)
1086 {
1087 double temt1rho01 = t1rho01[l]*G[l][l]*E[i][j][k][l];
1088 double temt1rho02 = t1rho02[l]*G[l][l]*E[i][j][k][l];
1089 double temt1rhop1 = t1rhop1[l]*G[l][l]*E[i][j][k][l];
1090 double temt1rhop2 = t1rhop2[l]*G[l][l]*E[i][j][k][l];
1091 double temt1rhom1 = t1rhom1[l]*G[l][l]*E[i][j][k][l];
1092 double temt1rhom2 = t1rhom2[l]*G[l][l]*E[i][j][k][l];
1093
1094 for(int m=0; m<4; m++)
1095 {
1096 double temt1kst1 = t1kst1[m]*G[m][m];
1097 double temt1kst2 = t1kst2[m]*G[m][m];
1098
1099 for(int n=0; n<4; n++)
1100 {
1101 if((i==j)||(i==k)||(i==l)||(j==k)||(j==l)||(k==l)) continue;
1102
1103 temp_PDF11 += t2D0_omega2kst1[m][n]*temt1kst1*p1omega02[n][i]*tempomega2*temt1omega02*temt1rho01*G[n][n];
1104 temp_PDF12 += t2D0_omega2kst1[m][n]*temt1kst1*p1omegap2[n][i]*tempomega2*temt1omegap2*temt1rhop2*G[n][n];
1105 temp_PDF13 += t2D0_omega2kst1[m][n]*temt1kst1*p1omegam2[n][i]*tempomega2*temt1omegam2*temt1rhom2*G[n][n];
1106 temp_PDF21 += t2D0_omega1kst2[m][n]*temt1kst2*p1omega01[n][i]*tempomega1*temt1omega01*temt1rho02*G[n][n];
1107 temp_PDF22 += t2D0_omega1kst2[m][n]*temt1kst2*p1omegap1[n][i]*tempomega1*temt1omegap1*temt1rhop1*G[n][n];
1108 temp_PDF23 += t2D0_omega1kst2[m][n]*temt1kst2*p1omegam1[n][i]*tempomega1*temt1omegam1*temt1rhom1*G[n][n];
1109 }
1110 }
1111 }
1112 }
1113 }
1114 }
1115
1116 temp_PDF11 = temp_PDF11*B2_D0_omega2*B1_omega20*B1_kst1*B1_V02;
1117 temp_PDF12 = temp_PDF12*B2_D0_omega2*B1_omega2p*B1_kst1*B1_Vp2;
1118 temp_PDF13 = temp_PDF13*B2_D0_omega2*B1_omega2m*B1_kst1*B1_Vm2;
1119 temp_PDF21 = temp_PDF21*B2_D0_omega1*B1_omega10*B1_kst2*B1_V01;
1120 temp_PDF22 = temp_PDF22*B2_D0_omega1*B1_omega1p*B1_kst2*B1_Vp1;
1121 temp_PDF23 = temp_PDF23*B2_D0_omega1*B1_omega1m*B1_kst2*B1_Vm1;
1122 }
1123
1124 //omega
1125 if(g0[i]==1){
1126 propagatorRBW(mass1sq,mass1[i],width1[i],somega1,srho01,sPi01,rRes2,1,pro0_01);//ome1--rho0 pi01
1127 propagatorRBW(mass1sq,mass1[i],width1[i],somega1,srhop1,sPim,rRes2,1,pro0_p1);//ome1--rhop (pi01)
1128 propagatorRBW(mass1sq,mass1[i],width1[i],somega1,srhom1,sPip,rRes2,1,pro0_m1);//ome1--rhom (pi01)
1129 propagatorRBW(mass1sq,mass1[i],width1[i],somega2,srho02,sPi02,rRes2,1,pro0_02);//ome2--rho0 pi02
1130 propagatorRBW(mass1sq,mass1[i],width1[i],somega2,srhop2,sPim,rRes2,1,pro0_p2);//ome2--rhop (pi02)
1131 propagatorRBW(mass1sq,mass1[i],width1[i],somega2,srhom2,sPip,rRes2,1,pro0_m2);//ome2--rhom(pi02)
1132 /* propagator(mass1sq,mass1[i],width1[i],somega1,pro0_01);//ome1--rho0 pi01
1133 propagator(mass1sq,mass1[i],width1[i],somega1,pro0_p1);//ome1--rhop (pi01)
1134 propagator(mass1sq,mass1[i],width1[i],somega1,pro0_m1);//ome1--rhom (pi01)
1135 propagator(mass1sq,mass1[i],width1[i],somega2,pro0_02);//ome2--rho0 pi02
1136 propagator(mass1sq,mass1[i],width1[i],somega2,pro0_p2);//ome2--rhop (pi02)
1137 propagator(mass1sq,mass1[i],width1[i],somega2,pro0_m2);
1138*/
1139 }
1140 else if(g0[i]==0){
1141 pro0_01[0] = 1; pro0_01[1] = 0; pro0_p1[0] = 1; pro0_p1[1] = 0; pro0_m1[0] = 1; pro0_m1[1] = 0;
1142 pro0_02[0] = 1; pro0_02[1] = 0; pro0_p2[0] = 1; pro0_p2[1] = 0; pro0_m2[0] = 1; pro0_m2[1] = 0;
1143 }
1144 //kst
1145 if(g1[i]==1){
1146 propagatorRBW(mass2sq,mass2[i],width2[i],skst1,sKs,sPi01,rRes2,1,pro1_1);//kst1
1147 propagatorRBW(mass2sq,mass2[i],width2[i],skst2,sKs,sPi02,rRes2,1,pro1_2);//kst2
1148 }
1149 else if(g1[i]==0){
1150 pro1_1[0] = 1; pro1_1[1] = 0;
1151 pro1_2[0] = 1; pro1_2[1] = 0;
1152 }
1153 //rho
1154 /* if(r1[i]==1){
1155 propagatorGS(0.60102807,0.77526,0.1478,srho01,sPip,sPim,rRes2,pro1V01);
1156 propagatorGS(0.60102807,0.77526,0.1478,srho02,sPip,sPim,rRes2,pro1V02);
1157 propagatorGS(0.60102807,0.77526,0.1478,srhop1,sPip,sPi01,rRes2,pro1Vp1);
1158 propagatorGS(0.60102807,0.77526,0.1478,srhop2,sPip,sPi02,rRes2,pro1Vp2);
1159 propagatorGS(0.60102807,0.77526,0.1478,srhom1,sPi01,sPim,rRes2,pro1Vm1);
1160 propagatorGS(0.60102807,0.77526,0.1478,srhom2,sPi02,sPim,rRes2,pro1Vm2);
1161
1162 }
1163 else if(r1[i]==0){
1164 pro1V01[0] = 1; pro1V01[1] = 0;
1165 pro1V02[0] = 1; pro1V02[1] = 0;
1166 pro1Vp1[0] = 1; pro1Vp1[1] = 0;
1167 pro1Vp2[0] = 1; pro1Vp2[1] = 0;
1168 pro1Vm1[0] = 1; pro1Vm1[1] = 0;
1169 pro1Vm2[0] = 1; pro1Vm2[1] = 0;
1170 }*/
1171
1172 Com_Multi(pro0_01,pro1V01,tpro_01); Com_Multi(tpro_01,pro1_2,pro_21);
1173 Com_Multi(pro0_p1,pro1Vp1,tpro_p1); Com_Multi(tpro_p1,pro1_2,pro_22);
1174 Com_Multi(pro0_m1,pro1Vm1,tpro_m1); Com_Multi(tpro_m1,pro1_2,pro_23);
1175
1176 Com_Multi(pro0_02,pro1V02,tpro_02); Com_Multi(tpro_02,pro1_1,pro_11);
1177 Com_Multi(pro0_p2,pro1Vp2,tpro_p2); Com_Multi(tpro_p2,pro1_1,pro_12);
1178 Com_Multi(pro0_m2,pro1Vm2,tpro_m2); Com_Multi(tpro_m2,pro1_1,pro_13);
1179
1180 amp_tmp1[0] =-1*(temp_PDF11*pro_11[0]) + temp_PDF12*pro_12[0] + temp_PDF13*pro_13[0];
1181 amp_tmp1[1] =-1*(temp_PDF11*pro_11[1]) + temp_PDF12*pro_12[1] + temp_PDF13*pro_13[1];
1182 amp_tmp2[0] =-1*(temp_PDF21*pro_21[0]) + temp_PDF22*pro_22[0] + temp_PDF23*pro_23[0];
1183 amp_tmp2[1] =-1*(temp_PDF21*pro_21[1]) + temp_PDF22*pro_22[1] + temp_PDF23*pro_23[1];
1184 // printf("shuchu1: %.20f,%lf,%.20f, %.20f\n", amp_tmp1[1],amp_tmp2[1],amp_tmp1[0],amp_tmp2[0]);
1185 }
1186 else if(modetype[i]== 40){
1187 if (g2[i]==0)//spin=0 D to rho omega
1188 {
1189 if(B1_omega10==-1) {
1190
1191 B1_omega10=BarrierN(1,somega1,srho01,sPi01,rRes2,mass1[i]);
1192 B1_omega1p=BarrierN(1,somega1,srhop1,sPim,rRes2,mass1[i]);
1193 B1_omega1m=BarrierN(1,somega1,srhom1,sPip,rRes2,mass1[i]);
1194
1195 B1_omega20=BarrierN(1,somega2,srho02,sPi02,rRes2,mass1[i]);
1196 B1_omega2p=BarrierN(1,somega2,srhop2,sPim,rRes2,mass1[i]);
1197 B1_omega2m=BarrierN(1,somega2,srhom2,sPip,rRes2,mass1[i]);
1198
1199 }
1200
1201 if(B1_D0_omega1==-1) {
1202 B1_D0_omega1=BarrierN(1,sD0,s1omega1,skpi02SW,rDs2,mD0M);
1203 B1_D0_omega2=BarrierN(1,sD0,s1omega2,skpi01SW,rDs2,mD0M);
1204 }
1205
1206 for(int i=0; i<4; i++)
1207 {
1208 for(int j=0; j<4; j++)
1209 {
1210 double tempomega1 = pomega1[j]*G[j][j]*G[i][i];// p(omega)
1211 double tempomega2 = pomega2[j]*G[j][j]*G[i][i];
1212
1213 for(int k=0; k<4; k++)
1214 {
1215 double temt1omega01 = t1omega01[k]*G[k][k];//Sn(ome)
1216 double temt1omega02 = t1omega02[k]*G[k][k];
1217 double temt1omegap1 = t1omegap1[k]*G[k][k];
1218 double temt1omegap2 = t1omegap2[k]*G[k][k];
1219 double temt1omegam1 = t1omegam1[k]*G[k][k];
1220 double temt1omegam2 = t1omegam2[k]*G[k][k];
1221
1222 for(int l=0; l<4; l++)
1223 {
1224 if((i==j)||(i==k)||(i==l)||(j==k)||(j==l)||(k==l)) continue;
1225
1226 double temt1rho01 = t1rho01[l]*G[l][l]*E[i][j][k][l];
1227 double temt1rho02 = t1rho02[l]*G[l][l]*E[i][j][k][l];
1228 double temt1rhop1 = t1rhop1[l]*G[l][l]*E[i][j][k][l];
1229 double temt1rhop2 = t1rhop2[l]*G[l][l]*E[i][j][k][l];
1230 double temt1rhom1 = t1rhom1[l]*G[l][l]*E[i][j][k][l];
1231 double temt1rhom2 = t1rhom2[l]*G[l][l]*E[i][j][k][l];
1232
1233 for(int m=0; m<4; m++)
1234 {
1235 temp_PDF11 += G[m][m]*t1D0KPi01SW[m]*p1omega02[m][i]*tempomega2*temt1omega02*temt1rho02;
1236 temp_PDF12 += G[m][m]*t1D0KPi01SW[m]*p1omegap2[m][i]*tempomega2*temt1omegap2*temt1rhop2;
1237 temp_PDF13 += G[m][m]*t1D0KPi01SW[m]*p1omegam2[m][i]*tempomega2*temt1omegam2*temt1rhom2;
1238
1239 temp_PDF21 += G[m][m]*t1D0KPi02SW[m]*p1omega01[m][i]*tempomega1*temt1omega01*temt1rho01;
1240 temp_PDF22 += G[m][m]*t1D0KPi02SW[m]*p1omegap1[m][i]*tempomega1*temt1omegap1*temt1rhop1;
1241 temp_PDF23 += G[m][m]*t1D0KPi02SW[m]*p1omegam1[m][i]*tempomega1*temt1omegam1*temt1rhom1;
1242 }
1243 }
1244 }
1245 }
1246 }
1247
1248 temp_PDF11 = temp_PDF11*B1_omega20*B1_D0_omega2*B1_V02;
1249 temp_PDF12 = temp_PDF12*B1_omega2p*B1_D0_omega2*B1_Vp2;
1250 temp_PDF13 = temp_PDF13*B1_omega2m*B1_D0_omega2*B1_Vm2;
1251 temp_PDF21 = temp_PDF21*B1_omega10*B1_D0_omega1*B1_V01;
1252 temp_PDF22 = temp_PDF22*B1_omega1p*B1_D0_omega1*B1_Vp1;
1253 temp_PDF23 = temp_PDF23*B1_omega1m*B1_D0_omega1*B1_Vm1;
1254 }
1255
1256 if(g0[i]==1){
1257 propagatorRBW(mass1sq,mass1[i],width1[i],somega1,srho01,sPi01,rRes2,1,pro0_01);//ome1--rho0 pi01
1258 propagatorRBW(mass1sq,mass1[i],width1[i],somega1,srhop1,sPim,rRes2,1,pro0_p1);//ome1--rhop (pi01)
1259 propagatorRBW(mass1sq,mass1[i],width1[i],somega1,srhom1,sPip,rRes2,1,pro0_m1);//ome1--rhom (pi01)
1260 propagatorRBW(mass1sq,mass1[i],width1[i],somega2,srho02,sPi02,rRes2,1,pro0_02);//ome2--rho0 pi02
1261 propagatorRBW(mass1sq,mass1[i],width1[i],somega2,srhop2,sPim,rRes2,1,pro0_p2);//ome2--rhop (pi02)
1262 propagatorRBW(mass1sq,mass1[i],width1[i],somega2,srhom2,sPip,rRes2,1,pro0_m2);//ome2--rhom(pi02)
1263 /* propagator(mass1sq,mass1[i],width1[i],somega1,pro0_01);//ome1--rho0 pi01
1264 propagator(mass1sq,mass1[i],width1[i],somega1,pro0_p1);//ome1--rhop (pi01)
1265 propagator(mass1sq,mass1[i],width1[i],somega1,pro0_m1);//ome1--rhom (pi01)
1266 propagator(mass1sq,mass1[i],width1[i],somega2,pro0_02);//ome2--rho0 pi02
1267 propagator(mass1sq,mass1[i],width1[i],somega2,pro0_p2);//ome2--rhop (pi02)
1268 propagator(mass1sq,mass1[i],width1[i],somega2,pro0_m2);
1269*/
1270 }
1271 else if(g0[i]==0){
1272 pro0_01[0] = 1; pro0_01[1] = 0; pro0_p1[0] = 1; pro0_p1[1] = 0; pro0_m1[0] = 1; pro0_m1[1] = 0;
1273 pro0_02[0] = 1; pro0_02[1] = 0; pro0_p2[0] = 1; pro0_p2[1] = 0; pro0_m2[0] = 1; pro0_m2[1] = 0;
1274 }
1275 //kst
1276 if(g1[i]==1){
1277 KPiSLASS(skpi01SW,sKs,sPi01,pro1_1);
1278 KPiSLASS(skpi02SW,sKs,sPi02,pro1_2);
1279 }
1280 else if(g1[i]==0){
1281 pro1_1[0] = 1; pro1_1[1] = 0;
1282 pro1_2[0] = 1; pro1_2[1] = 0;
1283 }
1284 //rho
1285 /* if(r1[i]==1){
1286 propagatorGS(0.60102807,0.77526,0.1478,srho01,sPip,sPim,rRes2,pro1V01);
1287 propagatorGS(0.60102807,0.77526,0.1478,srho02,sPip,sPim,rRes2,pro1V02);
1288 propagatorGS(0.60102807,0.77526,0.1478,srhop1,sPip,sPi01,rRes2,pro1Vp1);
1289 propagatorGS(0.60102807,0.77526,0.1478,srhop2,sPip,sPi02,rRes2,pro1Vp2);
1290 propagatorGS(0.60102807,0.77526,0.1478,srhom1,sPi01,sPim,rRes2,pro1Vm1);
1291 propagatorGS(0.60102807,0.77526,0.1478,srhom2,sPi02,sPim,rRes2,pro1Vm2);
1292
1293 }
1294 else if(r1[i]==0){
1295 pro1V01[0] = 1; pro1V01[1] = 0;
1296 pro1V02[0] = 1; pro1V02[1] = 0;
1297 pro1Vp1[0] = 1; pro1Vp1[1] = 0;
1298 pro1Vp2[0] = 1; pro1Vp2[1] = 0;
1299 pro1Vm1[0] = 1; pro1Vm1[1] = 0;
1300 pro1Vm2[0] = 1; pro1Vm2[1] = 0;
1301 }
1302 */
1303 Com_Multi(pro0_01,pro1V01,tpro_01); Com_Multi(tpro_01,pro1_2,pro_21);
1304 Com_Multi(pro0_p1,pro1Vp1,tpro_p1); Com_Multi(tpro_p1,pro1_2,pro_22);
1305 Com_Multi(pro0_m1,pro1Vm1,tpro_m1); Com_Multi(tpro_m1,pro1_2,pro_23);
1306
1307 Com_Multi(pro0_02,pro1V02,tpro_02); Com_Multi(tpro_02,pro1_1,pro_11);
1308 Com_Multi(pro0_p2,pro1Vp2,tpro_p2); Com_Multi(tpro_p2,pro1_1,pro_12);
1309 Com_Multi(pro0_m2,pro1Vm2,tpro_m2); Com_Multi(tpro_m2,pro1_1,pro_13);
1310
1311 amp_tmp1[0] =-1*(temp_PDF11*pro_11[0]) + temp_PDF12*pro_12[0] + temp_PDF13*pro_13[0];
1312 amp_tmp1[1] =-1*(temp_PDF11*pro_11[1]) + temp_PDF12*pro_12[1] + temp_PDF13*pro_13[1];
1313 amp_tmp2[0] =-1*(temp_PDF21*pro_21[0]) + temp_PDF22*pro_22[0] + temp_PDF23*pro_23[0];
1314 amp_tmp2[1] =-1*(temp_PDF21*pro_21[1]) + temp_PDF22*pro_22[1] + temp_PDF23*pro_23[1];
1315 }
1316
1317 else if (modetype[i]==4){
1318 if(g2[i]==0)
1319 {
1320 if(B1_omega10==-1) {
1321 B1_omega10=BarrierN(1,somega1,srho01,sPi01,rRes2,mass2[i]);
1322 B1_omega1p=BarrierN(1,somega1,srhop1,sPim,rRes2,mass2[i]);
1323 B1_omega1m=BarrierN(1,somega1,srhom1,sPip,rRes2,mass2[i]);
1324
1325 B1_omega20=BarrierN(1,somega2,srho02,sPi02,rRes2,mass2[i]);
1326 B1_omega2p=BarrierN(1,somega2,srhop2,sPim,rRes2,mass2[i]);
1327 B1_omega2m=BarrierN(1,somega2,srhom2,sPip,rRes2,mass2[i]);
1328
1329 }
1330
1331 if(B1_D0_b1==-1) {
1332 B1_D0_b1=BarrierN(1,sD0,sb10,sKs,rDs2,mD0M);
1333 }
1334
1335 for(int i=0; i<4; i++)
1336 {
1337 for(int j=0; j<4; j++)
1338 {
1339 double tempomega1 = pomega1[j]*G[j][j]*G[i][i];
1340 double tempomega2 = pomega2[j]*G[j][j]*G[i][i];
1341
1342 for(int k=0; k<4; k++)
1343 {
1344 double temt1omega01 = t1omega01[k]*G[k][k];
1345 double temt1omega02 = t1omega02[k]*G[k][k];
1346 double temt1omegap1 = t1omegap1[k]*G[k][k];
1347 double temt1omegap2 = t1omegap2[k]*G[k][k];
1348 double temt1omegam1 = t1omegam1[k]*G[k][k];
1349 double temt1omegam2 = t1omegam2[k]*G[k][k];
1350
1351 for(int l=0; l<4; l++)
1352 {
1353 double temt1rho01 = t1rho01[l]*G[l][l]*E[i][j][k][l];
1354 double temt1rho02 = t1rho02[l]*G[l][l]*E[i][j][k][l];
1355 double temt1rhom2 = t1rhom2[l]*G[l][l]*E[i][j][k][l];
1356 double temt1rhop1 = t1rhop1[l]*G[l][l]*E[i][j][k][l];
1357 double temt1rhop2 = t1rhop2[l]*G[l][l]*E[i][j][k][l];
1358 double temt1rhom1 = t1rhom1[l]*G[l][l]*E[i][j][k][l];
1359
1360 for(int m=0; m<4; m++)
1361 {
1362 double temt1D0_b1 = t1D0_b1[m]*G[m][m];
1363
1364 for(int n=0; n<4; n++)
1365 {
1366
1367 if((i==j)||(i==k)||(i==l)||(j==k)||(j==l)||(k==l)) continue;
1368 temp_PDF11 += temt1D0_b1*p1b1_2[m][n]*p1omega02[n][i]*tempomega2*temt1omega02*temt1rho02*G[n][n];
1369 temp_PDF12 += temt1D0_b1*p1b1_2[m][n]*p1omegap2[n][i]*tempomega2*temt1omegap2*temt1rhop2*G[n][n];
1370 temp_PDF13 += temt1D0_b1*p1b1_2[m][n]*p1omegam2[n][i]*tempomega2*temt1omegam2*temt1rhom2*G[n][n];
1371 temp_PDF21 += temt1D0_b1*p1b1_1[m][n]*p1omega01[n][i]*tempomega1*temt1omega01*temt1rho01*G[n][n];
1372 temp_PDF22 += temt1D0_b1*p1b1_1[m][n]*p1omegap1[n][i]*tempomega1*temt1omegap1*temt1rhop1*G[n][n];
1373 temp_PDF23 += temt1D0_b1*p1b1_1[m][n]*p1omegam1[n][i]*tempomega1*temt1omegam1*temt1rhom1*G[n][n];
1374 }
1375 }
1376 }
1377 }
1378 }
1379 }
1380 temp_PDF11 = temp_PDF11*B1_D0_b1*B1_omega20*B1_V02;
1381 temp_PDF12 = temp_PDF12*B1_D0_b1*B1_omega2p*B1_Vp2;
1382 temp_PDF13 = temp_PDF13*B1_D0_b1*B1_omega2m*B1_Vm2;
1383 temp_PDF21 = temp_PDF21*B1_D0_b1*B1_omega10*B1_V01;
1384 temp_PDF22 = temp_PDF22*B1_D0_b1*B1_omega1p*B1_Vp1;
1385 temp_PDF23 = temp_PDF23*B1_D0_b1*B1_omega1m*B1_Vm1;
1386 }
1387 else if(g2[i]==2)
1388 {
1389 if(B1_omega10==-1) {
1390
1391 B1_omega10=BarrierN(1,somega1,srho01,sPi01,rRes2,mass2[i]);
1392 B1_omega1p=BarrierN(1,somega1,srhop1,sPim,rRes2,mass2[i]);
1393 B1_omega1m=BarrierN(1,somega1,srhom1,sPip,rRes2,mass2[i]);
1394
1395 B1_omega20=BarrierN(1,somega2,srho02,sPi02,rRes2,mass2[i]);
1396 B1_omega2p=BarrierN(1,somega2,srhop2,sPim,rRes2,mass2[i]);
1397 B1_omega2m=BarrierN(1,somega2,srhom2,sPip,rRes2,mass2[i]);
1398 }
1399
1400 if(B1_D0_b1==-1) {
1401 B1_D0_b1=BarrierN(1,sD0,sb10,sKs,rDs2,mD0M);
1402 }
1403
1404 if(B2_b1_1==-1) {
1405 B2_b1_1=BarrierN(2,sb10,somega1,sPi02,rRes2,mass1[i]);
1406 B2_b1_2=BarrierN(2,sb10,somega2,sPi01,rRes2,mass1[i]);
1407 }
1408
1409 for(int i=0; i<4; i++)
1410 {
1411 for(int j=0; j<4; j++)
1412 {
1413 double tempomega1 = pomega1[j]*G[j][j]*G[i][i];
1414 double tempomega2 = pomega2[j]*G[j][j]*G[i][i];
1415
1416 for(int k=0; k<4; k++)
1417 {
1418 double temt1omega01 = t1omega01[k]*G[k][k];
1419 double temt1omega02 = t1omega02[k]*G[k][k];
1420 double temt1omegap1 = t1omegap1[k]*G[k][k];
1421 double temt1omegap2 = t1omegap2[k]*G[k][k];
1422 double temt1omegam1 = t1omegam1[k]*G[k][k];
1423 double temt1omegam2 = t1omegam2[k]*G[k][k];
1424
1425 for(int l=0; l<4; l++)
1426 {
1427 double temt1rho01 = t1rho01[l]*G[l][l]*E[i][j][k][l];
1428 double temt1rho02 = t1rho02[l]*G[l][l]*E[i][j][k][l];
1429 double temt1rhop1 = t1rhop1[l]*G[l][l]*E[i][j][k][l];
1430 double temt1rhop2 = t1rhop2[l]*G[l][l]*E[i][j][k][l];
1431 double temt1rhom1 = t1rhom1[l]*G[l][l]*E[i][j][k][l];
1432 double temt1rhom2 = t1rhom2[l]*G[l][l]*E[i][j][k][l];
1433
1434 for(int m=0; m<4; m++)
1435 {
1436 double temt1D0_b1 = t1D0_b1[m]*G[m][m];
1437
1438 for(int n=0; n<4; n++)
1439 {
1440 if((i==j)||(i==k)||(i==l)||(j==k)||(j==l)||(k==l)) continue;
1441 temp_PDF11 += temt1D0_b1*t2b1_2[m][n]*p1omega02[n][i]*tempomega2*temt1omega02*temt1rho02*G[n][n];
1442 temp_PDF12 += temt1D0_b1*t2b1_2[m][n]*p1omegap2[n][i]*tempomega2*temt1omegap2*temt1rhop2*G[n][n];
1443 temp_PDF13 += temt1D0_b1*t2b1_2[m][n]*p1omegam2[n][i]*tempomega2*temt1omegam2*temt1rhom2*G[n][n];
1444 temp_PDF21 += temt1D0_b1*t2b1_1[m][n]*p1omega01[n][i]*tempomega1*temt1omega01*temt1rho01*G[n][n];
1445 temp_PDF22 += temt1D0_b1*t2b1_1[m][n]*p1omegap1[n][i]*tempomega1*temt1omegap1*temt1rhop1*G[n][n];
1446 temp_PDF23 += temt1D0_b1*t2b1_1[m][n]*p1omegam1[n][i]*tempomega1*temt1omegam1*temt1rhom1*G[n][n];
1447 }
1448 }
1449 }
1450 }
1451 }
1452 }
1453 temp_PDF11 = temp_PDF11*B1_D0_b1*B2_b1_2*B1_omega20*B1_V02;
1454 temp_PDF12 = temp_PDF12*B1_D0_b1*B2_b1_2*B1_omega2p*B1_Vp2;
1455 temp_PDF13 = temp_PDF13*B1_D0_b1*B2_b1_2*B1_omega2m*B1_Vm2;
1456
1457 temp_PDF21 = temp_PDF21*B1_D0_b1*B2_b1_1*B1_omega10*B1_V01;
1458 temp_PDF22 = temp_PDF22*B1_D0_b1*B2_b1_1*B1_omega1p*B1_Vp1;
1459 temp_PDF23 = temp_PDF23*B1_D0_b1*B2_b1_1*B1_omega1m*B1_Vm1;
1460 }
1461 if(g0[i]==1) { //b1 omega
1462
1463 propagatorRBW(mass1sq,mass1[i],width1[i],sb10,somega1,sPi02,rRes2,g2[i],prob1_1);
1464 propagatorRBW(mass1sq,mass1[i],width1[i],sb10,somega2,sPi01,rRes2,g2[i],prob1_2);
1465
1466 } else if(g0[i]==0) {
1467 prob1_1[0] = 1; prob1_1[1] = 0;
1468 prob1_2[0] = 1; prob1_2[1] = 0;
1469 }
1470 if(g1[i]==1) { // omega--rho
1471
1472 propagatorRBW(mass2sq,mass2[i],width2[i],somega1,srho01,sPi01,rRes2,1,pro1_01);//01
1473 propagatorRBW(mass2sq,mass2[i],width2[i],somega1,srhop1,sPim,rRes2,1,pro1_p1);//p1
1474 propagatorRBW(mass2sq,mass2[i],width2[i],somega1,srhom1,sPip,rRes2,1,pro1_m1);//m1
1475 propagatorRBW(mass2sq,mass2[i],width2[i],somega2,srho02,sPi02,rRes2,1,pro1_02);//02
1476 propagatorRBW(mass2sq,mass2[i],width2[i],somega2,srhop2,sPim,rRes2,1,pro1_p2);//p2
1477 propagatorRBW(mass2sq,mass2[i],width2[i],somega2,srhom2,sPip,rRes2,1,pro1_m2);//m2
1478
1479 /* propagator(mass1sq,mass1[i],width1[i],somega1,pro0_01);//ome1--rho0 pi01
1480 propagator(mass1sq,mass1[i],width1[i],somega1,pro0_p1);//ome1--rhop (pi01)
1481 propagator(mass1sq,mass1[i],width1[i],somega1,pro0_m1);//ome1--rhom (pi01)
1482 propagator(mass1sq,mass1[i],width1[i],somega2,pro0_02);//ome2--rho0 pi02
1483 propagator(mass1sq,mass1[i],width1[i],somega2,pro0_p2);//ome2--rhop (pi02)
1484 propagator(mass1sq,mass1[i],width1[i],somega2,pro0_m2);
1485*/
1486 } else if(g1[i]==0) {
1487
1488 pro1_01[0] = 1; pro1_01[1] = 0; pro1_p1[0] = 1; pro1_p1[1] = 0; pro1_m1[0] = 1; pro1_m1[1] = 0;
1489 pro1_02[0] = 1; pro1_02[1] = 0; pro1_p2[0] = 1; pro1_p2[1] = 0; pro1_m2[0] = 1; pro1_m2[1] = 0;
1490
1491 }
1492
1493 //rho--pipi
1494 /* if(r1[i]==1){
1495 propagatorGS(0.60102807,0.77526,0.1478,srho01,sPip,sPim,rRes2,pro1V01);
1496 propagatorGS(0.60102807,0.77526,0.1478,srho02,sPip,sPim,rRes2,pro1V02);
1497 propagatorGS(0.60102807,0.77526,0.1478,srhop1,sPip,sPi01,rRes2,pro1Vp1);
1498 propagatorGS(0.60102807,0.77526,0.1478,srhop2,sPip,sPi02,rRes2,pro1Vp2);
1499 propagatorGS(0.60102807,0.77526,0.1478,srhom1,sPi01,sPim,rRes2,pro1Vm1);
1500 propagatorGS(0.60102807,0.77526,0.1478,srhom2,sPi02,sPim,rRes2,pro1Vm2);
1501
1502 }
1503 else if(r1[i]==0){
1504 pro1V01[0] = 1; pro1V01[1] = 0;
1505 pro1V02[0] = 1; pro1V02[1] = 0;
1506 pro1Vp1[0] = 1; pro1Vp1[1] = 0;
1507 pro1Vp2[0] = 1; pro1Vp2[1] = 0;
1508 pro1Vm1[0] = 1; pro1Vm1[1] = 0;
1509 pro1Vm2[0] = 1; pro1Vm2[1] = 0;
1510 }
1511 */
1512 Com_Multi(prob1_1,pro1_01,tpro_01); Com_Multi(tpro_01,pro1V01,pro_21);//10
1513 Com_Multi(prob1_1,pro1_p1,tpro_p1); Com_Multi(tpro_p1,pro1Vp1,pro_22);//1p
1514 Com_Multi(prob1_1,pro1_m1,tpro_m1); Com_Multi(tpro_m1,pro1Vm1,pro_23);//1m
1515
1516 Com_Multi(prob1_2,pro1_02,tpro_02); Com_Multi(tpro_02,pro1V02,pro_11);//20
1517 Com_Multi(prob1_2,pro1_p2,tpro_p2); Com_Multi(tpro_p2,pro1Vp2,pro_12);//2p
1518 Com_Multi(prob1_2,pro1_m2,tpro_m2); Com_Multi(tpro_m2,pro1Vm2,pro_13);//2m
1519
1520 amp_tmp1[0] = -1*(temp_PDF11*pro_11[0]) + temp_PDF12*pro_12[0] + temp_PDF13*pro_13[0];//b1--omega2
1521 amp_tmp1[1] = -1*(temp_PDF11*pro_11[1]) + temp_PDF12*pro_12[1] + temp_PDF13*pro_13[1];
1522 amp_tmp2[0] = -1*(temp_PDF21*pro_21[0]) + temp_PDF22*pro_22[0] + temp_PDF23*pro_23[0];//b1--omega1
1523 amp_tmp2[1] = -1*(temp_PDF21*pro_21[1]) + temp_PDF22*pro_22[1] + temp_PDF23*pro_23[1];
1524 // printf("shuchu2: %.20f,%.20f,%.20f, %.20f\n", amp_tmp1[1],amp_tmp2[1],amp_tmp1[0],amp_tmp2[0]);
1525 }
1526 //K*- a1+ --> rho+ pi0
1527 else if(modetype[i]==18){
1528 if(g2[i]==1020){ //0101---SS
1529
1530 if(B1_kst==-1){
1531 B1_kstm1 = BarrierN(1,skstm1,sKs,sPim,rRes2,mass2[i]); //m2==kst Fn(kst)
1532 B1_kstm2 = BarrierN(1,skstm2,sKs,sPim,rRes2,mass2[i]);
1533 }
1534 //P(a1)
1535 calG2(prhop1,Pi02,p1a1p_1);
1536 calG2(prhop2,Pi01,p1a1p_2);
1537
1538 calt1(PKs,Pim,t1kstm1);
1539 calt1(PKs,Pim,t1kstm2);
1540
1541 //t1rhop1 t1rhop2
1542
1543 for(int i=0; i<4; i++)
1544 {
1545 double t1_a1p_rhop1 = t1rhop1[i]*G[i][i];
1546 double t1_a1p_rhop2 = t1rhop2[i]*G[i][i];
1547
1548 for(int j=0; j<4; j++)
1549 {
1550 temp_PDF1 += t1_a1p_rhop1*t1kstm1[j]*G[j][j]*p1a1p_1[i][j];
1551 temp_PDF2 += t1_a1p_rhop2*t1kstm2[j]*G[j][j]*p1a1p_2[i][j];
1552 }
1553
1554 }
1555
1556 temp_PDF1 = temp_PDF1*B1_kstm1*B1_Vp1;
1557 temp_PDF2 = temp_PDF2*B1_kstm2*B1_Vp2;
1558
1559 }
1560/////////////////////////////SD
1561 else if(g2[i]==1022){
1562 if(B1_kstm==-1){
1563 B1_kstm1 = BarrierN(1,skstm1,sKs,sPim,rRes2,mass2[i]); //m2==kst Fn(kst)
1564 B1_kstm2 = BarrierN(1,skstm2,sKs,sPim,rRes2,mass2[i]);
1565 }
1566 if(B2_a1p_rhop1==-1){
1567 B2_a1p_rhop1 = BarrierN(2,sa1p,srhop1,sPi02,rRes2,mass1[i]);
1568 B2_a1p_rhop2 = BarrierN(2,sa1p,srhop2,sPi01,rRes2,mass1[i]);
1569 }
1570
1571 calt1(PKs,Pim,t1kstm1);
1572 calt1(PKs,Pim,t1kstm2);
1573 calt2(prhop1,Pi02,t2_a1p1);
1574 calt2(prhop2,Pi01,t2_a1p2);
1575 //t1rhop1 t1rhop2
1576
1577
1578 for(int i=0; i<4; i++)
1579 {
1580 double t1_a1p_rhop1 = t1rhop1[i]*G[i][i];
1581 double t1_a1p_rhop2 = t1rhop2[i]*G[i][i];
1582
1583 for(int j=0; j<4; j++)
1584 {
1585 temp_PDF1 += t1_a1p_rhop1*t1kstm1[j]*G[j][j]*t2_a1p1[j][i];
1586 temp_PDF2 += t1_a1p_rhop2*t1kstm2[j]*G[j][j]*t2_a1p2[j][i];
1587 }
1588
1589 }
1590 temp_PDF1 = temp_PDF1*B1_kstm1*B1_Vp1*B2_a1p_rhop1;
1591 temp_PDF2 = temp_PDF2*B1_kstm2*B1_Vp2*B2_a1p_rhop2;
1592
1593 }
1594
1595 else if(g2[i]==1120){
1596 if(B1_D0_a1pkstm==-1){
1597 B1_D0_a1pkstm1 = BarrierN(1,sD0,sa1p,skstm1,rDs2,mD0M);
1598 B1_D0_a1pkstm2 = BarrierN(1,sD0,sa1p,skstm2,rDs2,mD0M);
1599 }
1600
1601 if(B1_kstm==-1){
1602 B1_kstm1 = BarrierN(1,skstm1,sKs,sPim,rRes2,mass2[i]); //m2==kst Fn(kst)
1603 B1_kstm2 = BarrierN(1,skstm2,sKs,sPim,rRes2,mass2[i]);
1604 }
1605 //B1_Vp1 B1_Vp2--rho+
1606
1607 //Sn
1608 calt1(PKs,Pim,t1kstm1);
1609 calt1(PKs,Pim,t1kstm2);
1610 //t1rhop1 t1rhop2
1611 calt1(pa1p,pkstm1,t1D0_a1p1);
1612 calt1(pa1p,pkstm2,t1D0_a1p2);
1613
1614 //P
1615 calG2(prhop1,Pi02,p1a1p_1);
1616 calG2(prhop2,Pi01,p1a1p_2);
1617
1618 for(int i=0; i<4; i++)
1619 {
1620 double t1_a1p_kstm1 = t1kstm1[i]*G[i][i];
1621 double t1_a1p_kstm2 = t1kstm2[i]*G[i][i];
1622 for(int j=0; j<4; j++)
1623 {
1624 double t1_a1p_rhop1 = t1rhop1[j]*G[j][j];
1625 double t1_a1p_rhop2 = t1rhop2[i]*G[j][j];
1626
1627 for(int k=0; k<4; k++)
1628 {
1629 double t1_D0_a1pkstm1 = t1D0_a1p1[k]*G[k][k];
1630 double t1_D0_a1pkstm2 = t1D0_a1p2[k]*G[k][k];
1631
1632 for(int l=0; l<4; l++)
1633 {
1634 double pD0_a1pkstm1 = pD0[l]*G[l][l];
1635 double pD0_a1pkstm2 = pD0[l]*G[l][l];
1636
1637 for(int m=0; m<4; m++)
1638 {
1639 temp_PDF1 += t1_a1p_rhop1*t1_a1p_kstm1*t1_D0_a1pkstm1*pD0_a1pkstm1*p1a1p_1[m][j]*G[m][m]*E[i][k][l][m];
1640 temp_PDF2 += t1_a1p_rhop2*t1_a1p_kstm2*t1_D0_a1pkstm2*pD0_a1pkstm2*p1a1p_2[m][j]*G[m][m]*E[i][k][l][m];
1641 }
1642 }
1643 }
1644 }
1645 }
1646 temp_PDF1 = temp_PDF1*B1_D0_a1pkstm1*B1_kstm1*B1_Vp1;
1647 temp_PDF2 = temp_PDF2*B1_D0_a1pkstm2*B1_kstm2*B1_Vp2;
1648 }
1649
1650 else if(g2[i]==1122){//1121
1651 if(B1_D0_a1pkstm==-1){
1652 B1_D0_a1pkstm = BarrierN(1,sD0,sa1p,skstm,rDs2,mD0M);
1653 }
1654
1655 if(B1_kstm==-1){
1656 B1_kstm = BarrierN(1,skstm,sKs,sPim,rRes2,mass2[i]); //m2==kst Fn(kst)
1657 }
1658 //B1_Vp1 B1_Vp2--rho+
1659 if(B2_a1p_rhop1==-1){
1660 B2_a1p_rhop1 = BarrierN(2,sa1p,srhop1,sPi02,rRes2,mass1[i]);
1661 B2_a1p_rhop2 = BarrierN(2,sa1p,srhop2,sPi01,rRes2,mass1[i]);
1662 }
1663
1664 calt1(PKs,Pim,t1kstm);
1665 calt1(pa1p,pkstm,t1D0_a1p);
1666 //t1rhop1 t1rhop2
1667 calt2(prhop1,Pi02,t2_a1p1);
1668 calt2(prhop2,Pi01,t2_a1p2);
1669
1670
1671 for(int i=0; i<4; i++)
1672 {
1673 double t1_a1p_rhop1 = t1rhop1[i]*G[i][i];
1674 double t1_a1p_rhop2 = t1rhop2[i]*G[i][i];
1675
1676 for(int j=0; j<4; j++)
1677 {
1678 double t1_a1p_kstm = t1kstm[j]*G[j][j];
1679
1680 for(int k=0; k<4; k++)
1681 {
1682 double t1_D0_a1pkstm = t1D0_a1p[k]*G[k][k];
1683
1684 for(int l=0; l<4; l++)
1685 {
1686 double pD0_a1pkstm = pD0[l]*G[l][l];
1687
1688 for(int m=0; m<4; m++)
1689 {
1690 temp_PDF1 += t1_a1p_rhop1*t1_a1p_kstm*t1_D0_a1pkstm*pD0_a1pkstm*t2_a1p1[m][j]*G[m][m]*E[i][k][l][m];
1691 temp_PDF2 += t1_a1p_rhop2*t1_a1p_kstm*t1_D0_a1pkstm*pD0_a1pkstm*t2_a1p2[m][j]*G[m][m]*E[i][k][l][m];
1692 }
1693 }
1694 }
1695 }
1696 }
1697
1698 temp_PDF1 = temp_PDF1*B1_D0_a1pkstm*B1_kstm*B1_Vp1*B2_a1p_rhop1;
1699 temp_PDF2 = temp_PDF2*B1_D0_a1pkstm*B1_kstm*B1_Vp2*B2_a1p_rhop2;
1700
1701 }
1702
1703
1704 else if(g2[i]==1220){//2101
1705 if(B2_D0_a1pkstm==-1){
1706 B2_D0_a1pkstm = BarrierN(2,sD0,sa1p,skstm,rDs2,mD0M);
1707 }
1708 if(B1_kstm==-1){
1709 B1_kstm = BarrierN(1,skstm,sKs,sPim,rRes2,mass2[i]); //m2==kst Fn(kst)
1710 }
1711 //B1_Vp1 B1_Vp2--rho+
1712
1713 calt1(PKs,Pim,t1kstm);
1714 //t1rhop1 t1rhop2
1715 calt2(pa1p,pkstm,t2D0_a1p);
1716
1717 calG2(prhop1,Pi02,p1a1p_1);
1718 calG2(prhop2,Pi01,p1a1p_2);
1719
1720
1721 for(int i=0; i<4;i++)
1722 {
1723 double t1_a1p_rhop1 = t1rhop1[i]*G[i][i];
1724 double t1_a1p_rhop2 = t1rhop2[i]*G[i][i];
1725
1726 for(int j=0; j<4; j++)
1727 {
1728 double t1_a1p_kstm = t1kstm[j]*G[j][j];
1729
1730 for(int k=0; k<4; k++)
1731 {
1732 temp_PDF1 += t1_a1p_rhop1*t1_a1p_kstm*t2D0_a1p[k][i]*p1a1p_1[k][j];
1733 temp_PDF2 += t1_a1p_rhop2*t1_a1p_kstm*t2D0_a1p[k][i]*p1a1p_2[k][j];
1734 }
1735 }
1736 }
1737
1738 temp_PDF1 = temp_PDF1*B2_D0_a1pkstm*B1_kstm*B1_Vp1;
1739 temp_PDF2 = temp_PDF2*B2_D0_a1pkstm*B1_kstm*B1_Vp2;
1740
1741 }
1742
1743 else if(g2[i]==1222){//2121
1744 if(B2_D0_a1pkstm==-1){
1745 B2_D0_a1pkstm = BarrierN(2,sD0,sa1p,skstm,rDs2,mD0M);
1746 }
1747 if(B1_kstm==-1){
1748 B1_kstm = BarrierN(1,skstm,sKs,sPim,rRes2,mass2[i]); //m2==kst Fn(kst)
1749 }
1750 if(B2_a1p_rhop1==-1){
1751 B2_a1p_rhop1 = BarrierN(2,sa1p,srhop1,sPi02,rRes2,mass1[i]);
1752 B2_a1p_rhop2 = BarrierN(2,sa1p,srhop2,sPi01,rRes2,mass1[i]);
1753 }
1754 //B1_Vp1 B1_Vp2--rho+
1755
1756
1757 calt2(pa1p,pkstm,t2D0_a1p);
1758 calt1(PKs,Pim,t1kstm);
1759 calt2(prhop1,Pi02,t2_a1p1);
1760 calt2(prhop2,Pi01,t2_a1p2);
1761 //t1rhop1 t1rhop2
1762
1763
1764 for(int i=0; i<4;i++)
1765 {
1766 double t1_a1p_rhop1 = t1rhop1[i]*G[i][i];
1767 double t1_a1p_rhop2 = t1rhop2[i]*G[i][i];
1768
1769 for(int j=0; j<4; j++)
1770 {
1771 double t1_a1p_kstm = t1kstm[j]*G[j][j];
1772
1773 for(int k=0; k<4; k++)
1774 {
1775 temp_PDF1 += t1_a1p_rhop1*t1_a1p_kstm*t2D0_a1p[k][i]*t2_a1p1[k][j]*G[k][k];
1776 temp_PDF2 += t1_a1p_rhop2*t1_a1p_kstm*t2D0_a1p[k][i]*t2_a1p2[k][j]*G[k][k];
1777 }
1778 }
1779 }
1780
1781 temp_PDF1 = temp_PDF1*B2_D0_a1pkstm*B1_kstm*B1_Vp1*B2_a1p_rhop1;
1782 temp_PDF2 = temp_PDF2*B2_D0_a1pkstm*B1_kstm*B1_Vp2*B2_a1p_rhop2;
1783
1784 }
1785
1786
1787 if(g0[i] == 1){//a1
1788 if(g2[i]==1020||g2[i]==1120||g2[i]==1220) {
1789 propagatorRBW(mass1sq,mass1[i],width1[i],sa1p,srhop1,sPi02,rRes2,0,pro0_p1);
1790 propagatorRBW(mass1sq,mass1[i],width1[i],sa1p,srhop2,sPi01,rRes2,0,pro0_p2);
1791 }
1792 else{
1793 propagatorRBW(mass1sq,mass1[i],width1[i],sa1p,srhop1,sPi02,rRes2,2,pro0_p1);
1794 propagatorRBW(mass1sq,mass1[i],width1[i],sa1p,srhop2,sPi01,rRes2,2,pro0_p2);
1795
1796 }
1797
1798 }
1799 else if(g0[i]==0){
1800 pro0_p1[0] = 1; pro0_p1[1] = 0;
1801 pro0_p2[0] = 1; pro0_p2[1] = 0;
1802 }
1803
1804
1805 //kst
1806 if(g1[i]==1){
1807 propagatorRBW(mass2sq,mass2[i],width2[i],skstm1,sKs,sPim,rRes2,1,pro1_1);//kst1
1808 propagatorRBW(mass2sq,mass2[i],width2[i],skstm2,sKs,sPim,rRes2,1,pro1_2);//kst1
1809 // KPiSLASS(skstm,sKs,sPim,pro1_1);
1810 // KPiSLASS(skstm,sKs,sPim,pro1_2);
1811 }
1812 else if(g1[i]==0){
1813 pro1_1[0] = 1; pro1_1[1] = 0;
1814 pro1_2[0] = 1; pro1_2[1] = 0;
1815 }
1816
1817 //rho
1818 /* if(r1[i]==1){
1819 propagatorGS(0.60102807,0.77526,0.1478,srhop1,sPip,sPi01,rRes2,pro1Vp1);
1820 propagatorGS(0.60102807,0.77526,0.1478,srhop2,sPip,sPi02,rRes2,pro1Vp2);
1821
1822 }
1823 else if(r1[i]==0){
1824 pro1Vp1[0] = 1; pro1Vp1[1] = 0;
1825 pro1Vp2[0] = 1; pro1Vp2[1] = 0;
1826 }
1827 */
1828
1829 Com_Multi(pro0_p1,pro1Vp1,tpro_p1); Com_Multi(tpro_p1,pro1_1,pro_11);
1830 Com_Multi(pro0_p2,pro1Vp2,tpro_p2); Com_Multi(tpro_p2,pro1_2,pro_12);
1831
1832 amp_tmp1[0] = temp_PDF1*pro_11[0];
1833 amp_tmp1[1] = temp_PDF1*pro_11[1];
1834 amp_tmp2[0] = temp_PDF2*pro_12[0];
1835 amp_tmp2[1] = temp_PDF2*pro_12[1];
1836 }
1837
1838 else if(modetype[i]==19){
1839 if(g2[i]==1020){ //0101---SS
1840 if(B1_D0_a1pkpim1==-1){
1841 B1_D0_a1pkpim1 = BarrierN(1,sD0,sa1p,skstm1,rDs2,mD0M);
1842 B1_D0_a1pkpim2 = BarrierN(1,sD0,sa1p,skstm1,rDs2,mD0M);
1843 }
1844
1845 //P(a1)
1846 calG2(prhop1,Pi02,p1a1p_1);
1847 calG2(prhop2,Pi01,p1a1p_2);
1848
1849 calt1(pa1p,pkstm1,t1D0_kpimSW1);
1850 calt1(pa1p,pkstm1,t1D0_kpimSW2);
1851
1852 //t1rhop1 t1rhop2
1853
1854 for(int i=0; i<4; i++)
1855 {
1856 double t1_a1p_rhop1 = t1rhop1[i]*G[i][i];
1857 double t1_a1p_rhop2 = t1rhop2[i]*G[i][i];
1858
1859 for(int j=0; j<4; j++)
1860 {
1861 temp_PDF1 += t1_a1p_rhop1*t1D0_kpimSW1[j]*G[j][j]*p1a1p_1[i][j];
1862 temp_PDF2 += t1_a1p_rhop2*t1D0_kpimSW2[j]*G[j][j]*p1a1p_2[i][j];
1863 }
1864
1865 }
1866
1867 temp_PDF1 = temp_PDF1*B1_D0_a1pkpim1*B1_Vp1;
1868 temp_PDF2 = temp_PDF2*B1_D0_a1pkpim2*B1_Vp2;
1869
1870 }
1871 if(g0[i] == 1){//a1
1872 propagatorRBW(mass1sq,mass1[i],width1[i],sa1p,srhop1,sPi02,rRes2,0,pro0_p1);
1873 propagatorRBW(mass1sq,mass1[i],width1[i],sa1p,srhop2,sPi01,rRes2,0,pro0_p2);
1874 }
1875 else if(g0[i]==0){
1876 pro0_p1[0] = 1; pro0_p1[1] = 0;
1877 pro0_p2[0] = 1; pro0_p2[1] = 0;
1878 }
1879
1880 //KPISW
1881 if(g1[i]==1){
1882 KPiSLASS(skstm1,sKs,sPim,pro1_1);
1883 KPiSLASS(skstm1,sKs,sPim,pro1_2);
1884 }
1885 else if(g1[i]==0){
1886 pro1_1[0] = 1; pro1_1[1] = 0;
1887 }
1888
1889 //rho
1890 if(r1[i]==1){
1891 propagatorGS(0.60102807,0.77526,0.1478,srhop1,sPip,sPi01,rRes2,pro1Vp1);
1892 propagatorGS(0.60102807,0.77526,0.1478,srhop2,sPip,sPi02,rRes2,pro1Vp2);
1893
1894 }
1895 else if(r1[i]==0){
1896 pro1Vp1[0] = 1; pro1Vp1[1] = 0;
1897 pro1Vp2[0] = 1; pro1Vp2[1] = 0;
1898 }
1899
1900
1901 Com_Multi(pro0_p1,pro1Vp1,tpro_p1); Com_Multi(tpro_p1,pro1_1,pro_11);
1902 Com_Multi(pro0_p2,pro1Vp2,tpro_p2); Com_Multi(tpro_p2,pro1_2,pro_12);
1903
1904 amp_tmp1[0] = temp_PDF1*pro_11[0];
1905 amp_tmp1[1] = temp_PDF1*pro_11[1];
1906 amp_tmp2[0] = temp_PDF2*pro_12[0];
1907 amp_tmp2[1] = temp_PDF2*pro_12[1];
1908 }
1909
1910
1911 else if(modetype[i]==20){//a1_0 K*0 a1_0 to rho_0 pi0 K*0 to KS0 pi0 rho_0 to pi+ pi-
1912 if(g2[i]==1020){ //0101---SS
1913 if(B1_a10_kst0==-1){
1914 B1_kst1=BarrierN(1,skst1,sKs,sPi01,rRes2,mass2[i]);
1915 B1_kst2=BarrierN(1,skst2,sKs,sPi02,rRes2,mass2[i]);
1916 }
1917 calG2(prhop1,Pim,p1_a102);
1918 calG2(prhop2,Pim,p2_a102);
1919 for(int i=0; i<4; i++)
1920 {
1921 double temt1rhop1 = t1rhop1[i]*G[i][i];
1922 double temt1rhop2 = t1rhop2[i]*G[i][i];
1923 for(int j=0; j<4; j++)
1924 {
1925 temp_PDF1 += temt1rhop2*t1kst1[j]*G[j][j]*p2_a102[i][j];
1926 temp_PDF2 += temt1rhop1*t1kst2[j]*G[j][j]*p1_a102[i][j];
1927 }
1928 }
1929
1930 temp_PDF1 = temp_PDF1*B1_kst1*B1_Vp2;
1931 temp_PDF2 = temp_PDF2*B1_kst2*B1_Vp1;
1932 }
1933/////////////////////////////SD
1934 else if(g2[i]==1022){
1935 if(B1_a10_kst0==-1){
1936 B1_kst1=BarrierN(1,skst1,sKs,sPi01,rRes2,mass2[i]);
1937 B1_kst2=BarrierN(1,skst2,sKs,sPi02,rRes2,mass2[i]);
1938 }
1939 if(B2_a10_rho==-1){
1940 B1_a11p=BarrierN(2,sa101,srhop1,sPim,rRes2,mass1[i]);
1941 B2_a12p=BarrierN(2,sa102,srhop2,sPim,rRes2,mass1[i]);
1942 }
1943
1944 calt2(prhop1,Pim,t1_a102);
1945 calt2(prhop2,Pim,t2_a102);
1946 //t1rhop1 t1rhop2
1947
1948 for(int i=0; i<4; i++)
1949 {
1950 double t1_a10_rhop1 = t1rhop1[i]*G[i][i];
1951 double t1_a10_rhop2 = t1rhop2[i]*G[i][i];
1952 for(int j=0; j<4; j++)
1953 {
1954 temp_PDF1 += t1_a10_rhop2*t1kst1[j]*G[j][j]*t2_a102[j][i];
1955 temp_PDF2 += t1_a10_rhop1*t1kst2[j]*G[j][j]*t1_a102[j][i];
1956 }
1957
1958 }
1959 temp_PDF1 = temp_PDF1*B1_kst1*B1_Vp2*B2_a12p;
1960 temp_PDF2 = temp_PDF2*B1_kst2*B1_Vp1*B1_a11p;
1961
1962 }
1963
1964 else if(g2[i]==1120){
1965 if(B1_D0_a10kst0==-1){
1966 B1_D0_a10kst01 = BarrierN(1,sD0,sa102,skst1,rDs2,mD0M);
1967 B1_D0_a10kst02 = BarrierN(1,sD0,sa101,skst2,rDs2,mD0M);
1968 }
1969 if(B1_a10_kst0==-1){
1970 B1_kst1=BarrierN(1,skst1,sKs,sPi01,rRes2,mass2[i]);
1971 B1_kst2=BarrierN(1,skst2,sKs,sPi02,rRes2,mass2[i]);
1972 }
1973 //B1_Vp1 B1_Vp2--rho+
1974
1975 //t1rhop1 t1rhop2
1976 calt1(pa101,pkst2,t1D0_a101);
1977 calt1(pa102,pkst1,t1D0_a102);
1978
1979 //P
1980 calG2(prhop1,Pim,p1a1p_1);
1981 calG2(prhop2,Pim,p1a1p_2);
1982 for(int i=0; i<4; i++)
1983 {
1984 double t1_a10_rhop1 = t1rhop1[i]*G[i][i];
1985 double t1_a10_rhop2 = t1rhop2[i]*G[i][i];
1986
1987 for(int j=0; j<4; j++)
1988 {
1989 double t1_a10_kst01 = t1kst1[j]*G[j][j];
1990 double t1_a10_kst02 = t1kst2[j]*G[j][j];
1991
1992 for(int k=0; k<4; k++)
1993 {
1994 double t1_D0_a10kst01 = t1D0_a101[k]*G[k][k];
1995 double t1_D0_a10kst02 = t1D0_a102[k]*G[k][k];
1996
1997 for(int l=0; l<4; l++)
1998 { double pD0_a10kst1 = pD0[l]*G[l][l];
1999 double pD0_a10kst2 = pD0[l]*G[l][l];
2000
2001 for(int m=0; m<4; m++)
2002 {
2003 temp_PDF1 += t1_a10_rhop2*t1_a10_kst01*t1_D0_a10kst02*pD0_a10kst1*p1a1p_2[m][j]*G[m][m]*E[i][k][l][m];
2004 temp_PDF2 += t1_a10_rhop1*t1_a10_kst02*t1_D0_a10kst01*pD0_a10kst2*p1a1p_1[m][j]*G[m][m]*E[i][k][l][m];
2005 }
2006 }
2007 }
2008 }
2009 }
2010 temp_PDF1 = temp_PDF1*B1_D0_a10kst01*B1_kst1*B1_Vp2;
2011 temp_PDF2 = temp_PDF2*B1_D0_a10kst02*B1_kst2*B1_Vp1;
2012
2013 }
2014
2015 else if(g2[i]==1122){//1121
2016 if(B1_D0_a10kst0==-1){
2017 B1_D0_a10kst01 = BarrierN(1,sD0,sa102,skst1,rDs2,mD0M);
2018 B1_D0_a10kst02 = BarrierN(1,sD0,sa101,skst2,rDs2,mD0M);
2019 }
2020
2021 if(B1_a10_kst0==-1){
2022 B1_kst1=BarrierN(1,skst1,sKs,sPi01,rRes2,mass2[i]);
2023 B1_kst2=BarrierN(1,skst2,sKs,sPi02,rRes2,mass2[i]);
2024 }
2025 //B1_Vp1 B1_Vp2--rho+
2026 if(B2_a10_rho==-1){
2027 B1_a11p=BarrierN(2,sa101,srhop1,sPim,rRes2,mass1[i]);
2028
2029 B2_a12p=BarrierN(2,sa102,srhop2,sPim,rRes2,mass1[i]);
2030 }
2031
2032 calt1(pa101,pkst2,t1D0_a101);
2033 calt1(pa102,pkst1,t1D0_a102);
2034 //t1rhop1 t1rhop2
2035 calt2(prhop1,Pim,t1_a102);
2036 calt2(prhop2,Pim,t2_a102);
2037 for(int i=0; i<4; i++)
2038 {
2039 double t1_a10_rhop1 = t1rhop1[i]*G[i][i];
2040 double t1_a10_rhop2 = t1rhop2[i]*G[i][i];
2041
2042 for(int j=0; j<4; j++)
2043 {
2044 double t1_a10_kst01 = t1kst1[j]*G[j][j];
2045 double t1_a10_kst02 = t1kst2[j]*G[j][j];
2046 for(int k=0; k<4; k++)
2047 {
2048 double t1_D0_a101kst02 = t1D0_a101[k]*G[k][k];
2049 double t1_D0_a102kst01 = t1D0_a102[k]*G[k][k];
2050
2051 for(int l=0; l<4; l++)
2052 {
2053 double pD0_a10kst0 = pD0[l]*G[l][l];
2054
2055 for(int m=0; m<4; m++)
2056 {
2057 temp_PDF1 += t1_a10_rhop2*t1_a10_kst01*t1_D0_a101kst02*pD0_a10kst0*t2_a102[m][j]*G[m][m]*E[i][k][l][m];
2058 temp_PDF2 += t1_a10_rhop1*t1_a10_kst02*t1_D0_a102kst01*pD0_a10kst0*t1_a102[m][j]*G[m][m]*E[i][k][l][m];
2059 }
2060 }
2061 }
2062
2063 }
2064 }
2065 temp_PDF1 = temp_PDF1*B1_D0_a10kst01*B1_kst1*B1_Vp2*B2_a12p;
2066
2067 temp_PDF2 = temp_PDF2*B1_D0_a10kst02*B1_kst2*B1_Vp1*B1_a11p;
2068
2069 }
2070
2071
2072 else if(g2[i]==1220){//2101
2073 if(B2_D0_a10kst0==-1){
2074 B2_D0_a10kst01 = BarrierN(2,sD0,sa102,skst1,rDs2,mD0M);
2075 B2_D0_a10kst02 = BarrierN(2,sD0,sa101,skst2,rDs2,mD0M);
2076 }
2077 if(B1_a10_kst0==-1){
2078 B1_kst1=BarrierN(1,skst1,sKs,sPi01,rRes2,mass2[i]);
2079 B1_kst2=BarrierN(1,skst2,sKs,sPi02,rRes2,mass2[i]);
2080 }
2081 //B1_Vp1 B1_Vp2--rho+
2082
2083 //t1rhop1 t1rhop2
2084 calt2(pa101,pkst2,t2D0_a101);
2085 calt2(pa102,pkst1,t2D0_a102);
2086
2087 calG2(prhop1,Pim,p1a1p_1);
2088 calG2(prhop2,Pim,p1a1p_2);
2089 for(int i=0; i<4;i++)
2090 {
2091 double t1_a10_rhop1 = t1rhop1[i]*G[i][i];
2092 double t1_a10_rhop2 = t1rhop2[i]*G[i][i];
2093
2094 for(int j=0; j<4; j++)
2095 {
2096 double t1_a10_kst01 = t1kst1[j]*G[j][j];
2097 double t1_a10_kst02 = t1kst2[j]*G[j][j];
2098
2099 for(int k=0; k<4; k++)
2100 {
2101 temp_PDF1 += t1_a10_rhop2*t1_a10_kst01*t2D0_a102[k][i]*p1a1p_2[k][j];
2102 temp_PDF2 += t1_a10_rhop1*t1_a10_kst02*t2D0_a101[k][i]*p1a1p_1[k][j];
2103 }
2104 }
2105 }
2106
2107 temp_PDF1 = temp_PDF1*B2_D0_a10kst01*B1_kst1*B1_Vp2;
2108 temp_PDF2 = temp_PDF2*B2_D0_a10kst02*B1_kst2*B1_Vp1;
2109
2110 }
2111
2112 else if(g2[i]==1222){//2121
2113 if(B2_D0_a10kst0==-1){
2114 B2_D0_a10kst01 = BarrierN(2,sD0,sa102,skst1,rDs2,mD0M);
2115 B2_D0_a10kst02 = BarrierN(2,sD0,sa101,skst2,rDs2,mD0M);
2116 }
2117 if(B1_a10_kst0==-1){
2118 B1_kst1=BarrierN(1,skst1,sKs,sPi01,rRes2,mass2[i]);
2119 B1_kst2=BarrierN(1,skst2,sKs,sPi02,rRes2,mass2[i]);
2120 }
2121 if(B2_a10_rho==-1){
2122 B1_a11p=BarrierN(2,sa101,srhop1,sPim,rRes2,mass1[i]);
2123 B2_a12p=BarrierN(2,sa102,srhop2,sPim,rRes2,mass1[i]);
2124 }
2125 //B1_Vp1 B1_Vp2--rho+
2126
2127
2128 calt2(pa101,pkst2,t2D0_a101);
2129 calt2(pa102,pkst1,t2D0_a102);
2130
2131 calt2(prhop1,Pim,t1_a102);
2132 calt2(prhop2,Pim,t2_a102);
2133 //t1rhop1 t1rhop2
2134
2135
2136 for(int i=0; i<4;i++)
2137 {
2138 double t1_a10_rhop1 = t1rhop1[i]*G[i][i];
2139 double t1_a10_rhop2 = t1rhop2[i]*G[i][i];
2140 for(int j=0; j<4; j++)
2141 {
2142 double t1_a10_kst01 = t1kst1[j]*G[j][j];
2143 double t1_a10_kst02 = t1kst2[j]*G[j][j];
2144
2145 for(int k=0; k<4; k++)
2146 {
2147 temp_PDF1 += t1_a10_rhop2*t1_a10_kst01*t2D0_a102[k][i]*t2_a102[k][j]*G[k][k];
2148 temp_PDF2 += t1_a10_rhop1*t1_a10_kst02*t2D0_a101[k][i]*t1_a102[k][j]*G[k][k];
2149 }
2150 }
2151 }
2152
2153 temp_PDF1 = temp_PDF1*B2_D0_a10kst01*B1_kst1*B1_Vp2*B2_a12p;
2154 temp_PDF2 = temp_PDF2*B2_D0_a10kst02*B1_kst2*B1_Vp1*B1_a11p;
2155
2156 }
2157
2158
2159 if(g0[i] == 1){//a1
2160 if(g2[i]==1020||g2[i]==1120||g2[i]==1220) {
2161 propagatorRBW(mass1sq,mass1[i],width1[i],sa101,srhop1,sPim,rRes2,0,proa_p1);
2162 propagatorRBW(mass1sq,mass1[i],width1[i],sa102,srhop2,sPim,rRes2,0,proa_p2);
2163 }
2164 else{
2165 propagatorRBW(mass1sq,mass1[i],width1[i],sa101,srhop1,sPim,rRes2,2,proa_p1);
2166 propagatorRBW(mass1sq,mass1[i],width1[i],sa102,srhop2,sPim,rRes2,2,proa_p2);
2167 }
2168 }
2169 else if(g0[i]==0){
2170 proa_p1[0] = 1; proa_p1[1] = 0;
2171 proa_p2[0] = 1; proa_p2[1] = 0;
2172 }
2173
2174
2175 if(g1[i]==1){
2176 propagatorRBW(mass2sq,mass2[i],width2[i],skst1,sKs,sPi01,rRes2,1,pro1_1);//kst1
2177 propagatorRBW(mass2sq,mass2[i],width2[i],skst2,sKs,sPi02,rRes2,1,pro1_2);//kst2
2178 }
2179 else if(g1[i]==0){
2180 pro1_1[0] = 1; pro1_1[1] = 0;
2181 pro1_2[0] = 1; pro1_2[1] = 0;
2182 }
2183
2184 /* if(r1[i]==1){
2185 propagatorGS(0.60102807,0.77526,0.1478,srhop1,sPip,sPi01,rRes2,pro1Vp1);
2186 propagatorGS(0.60102807,0.77526,0.1478,srhop2,sPip,sPi02,rRes2,pro1Vp2);
2187 }
2188 else if(r1[i]==0){
2189 pro1Vp1[0] = 1; pro1Vp1[1] = 0;
2190 pro1Vp2[0] = 1; pro1Vp2[1] = 0;
2191 }
2192 */
2193
2194 Com_Multi(proa_p1,pro1Vp1,tpro_p1); Com_Multi(tpro_p1,pro1_2,pro_12);//1p
2195 Com_Multi(proa_p2,pro1Vp2,tpro_p2); Com_Multi(tpro_p2,pro1_1,pro_22);//2p
2196
2197 amp_tmp1[0] = temp_PDF1*pro_22[0] ;//b1--omega2
2198 amp_tmp1[1] = temp_PDF1*pro_22[1] ;
2199 amp_tmp2[0] = temp_PDF2*pro_12[0] ;//b1--omega1
2200 amp_tmp2[1] = temp_PDF2*pro_12[1] ;
2201 }
2202 else if(modetype[i]==33){//a1_0 K*0 a1_0 to rho_0 pi0 K*0 to KS0 pi0 rho_0 to pi+ pi-
2203 if(g2[i]==1020){ //0101---SS
2204 if(B1_a10_kst==-1){
2205 B1_kst1=BarrierN(1,skst1,sKs,sPi01,rRes2,mass2[i]);
2206 B1_kst2=BarrierN(1,skst2,sKs,sPi02,rRes2,mass2[i]);
2207 }
2208 //P(a1)
2209 calG2(prhom1,Pip,p1_a103);
2210 calG2(prhom2,Pip,p2_a103);
2211 // printf("shuchu1: %lf, %lf\n", t1kst1[0],p2_a102[0][0]);
2212 for(int i=0; i<4; i++)
2213 {
2214 double temt1rhom1 = t1rhom1[i]*G[i][i];
2215 double temt1rhom2 = t1rhom2[i]*G[i][i];
2216 for(int j=0; j<4; j++)
2217 {
2218 temp_PDF1 += temt1rhom2*t1kst1[j]*G[j][j]*p2_a103[i][j];
2219 temp_PDF2 += temt1rhom1*t1kst2[j]*G[j][j]*p1_a103[i][j];
2220 // printf("shuchu1: %lf, %lf\n", p2_a103[0][0],p1_a103[0][0]);
2221 }
2222
2223 }
2224 // printf("shuchu1: %lf, %lf\n", B1_kst1,B1_Vm2);
2225 temp_PDF1 = temp_PDF1*B1_kst1*B1_Vm2;
2226 temp_PDF2 = temp_PDF2*B1_kst2*B1_Vm1;
2227 }
2228 else if(g2[i]==1022){
2229 if(B1_a10_kst==-1){
2230 B1_kst1=BarrierN(1,skst1,sKs,sPi01,rRes2,mass2[i]);
2231 B1_kst2=BarrierN(1,skst2,sKs,sPi02,rRes2,mass2[i]);
2232 }
2233 if(B2_a10_rho==-1){
2234 B1_a11m=BarrierN(2,sa101,srhom1,sPip,rRes2,mass1[i]);
2235
2236 B2_a12m=BarrierN(2,sa102,srhom2,sPip,rRes2,mass1[i]);
2237 }
2238
2239 calt2(prhom1,Pip,t1_a103);
2240 calt2(prhom2,Pip,t2_a103);
2241 //t1rhop1 t1rhop2
2242
2243 for(int i=0; i<4; i++)
2244 {
2245 double t1_a10_rhom1 = t1rhom1[i]*G[i][i];
2246 double t1_a10_rhom2 = t1rhom2[i]*G[i][i];
2247 for(int j=0; j<4; j++)
2248 {
2249 temp_PDF1 += t1_a10_rhom2*t1kst1[j]*G[j][j]*t2_a103[j][i];
2250
2251 temp_PDF2 += t1_a10_rhom1*t1kst2[j]*G[j][j]*t1_a103[j][i];
2252 }
2253
2254 }
2255 temp_PDF1 = temp_PDF1*B1_kst1*B1_Vm2*B2_a12m;
2256 temp_PDF2 = temp_PDF2*B1_kst2*B1_Vm1*B1_a11m;
2257
2258 }
2259 else if(g2[i]==1120){
2260 if(B1_D0_a10kst0==-1){
2261 B1_D0_a10kst01 = BarrierN(1,sD0,sa102,skst1,rDs2,mD0M);
2262 B1_D0_a10kst02 = BarrierN(1,sD0,sa101,skst2,rDs2,mD0M);
2263 }
2264 if(B1_a10_kst==-1){
2265 B1_kst1=BarrierN(1,skst1,sKs,sPi01,rRes2,mass2[i]);
2266 B1_kst2=BarrierN(1,skst2,sKs,sPi02,rRes2,mass2[i]);
2267
2268 }
2269 //B1_Vp1 B1_Vp2--rho+
2270
2271 //t1rhop1 t1rhop2
2272 calt1(pa101,pkst2,t1D0_a101);
2273 calt1(pa102,pkst1,t1D0_a102);
2274
2275 //P
2276 calG2(prhom1,Pip,p1a1m_1);
2277 calG2(prhom2,Pip,p1a1m_2);
2278 for(int i=0; i<4; i++)
2279 {
2280 double t1_a10_kst01 = t1kst1[i]*G[i][i];
2281 double t1_a10_kst02 = t1kst2[i]*G[i][i];
2282 for(int j=0; j<4; j++)
2283 {
2284 double t1_a10_rhom1 = t1rhom1[j]*G[j][j];
2285 double t1_a10_rhom2 = t1rhom2[j]*G[j][j];
2286 for(int k=0; k<4; k++)
2287 {
2288 double t1_D0_a10kst01 = t1D0_a101[k]*G[k][k];
2289 double t1_D0_a10kst02 = t1D0_a102[k]*G[k][k];
2290
2291 for(int l=0; l<4; l++)
2292 { double pD0_a10kst1 = pD0[l]*G[l][l];
2293 double pD0_a10kst2 = pD0[l]*G[l][l];
2294
2295 for(int m=0; m<4; m++)
2296 {
2297 temp_PDF1 += t1_a10_rhom2*t1_a10_kst01*t1_D0_a10kst02*pD0_a10kst1*p1a1m_2[m][j]*G[m][m]*E[i][k][l][m];
2298 temp_PDF2 += t1_a10_rhom1*t1_a10_kst02*t1_D0_a10kst01*pD0_a10kst2*p1a1m_1[m][j]*G[m][m]*E[i][k][l][m];
2299 }
2300 }
2301 }
2302 }
2303 }
2304 temp_PDF1 = temp_PDF1*B1_D0_a10kst01*B1_kst1*B1_Vm2;
2305 temp_PDF2 = temp_PDF2*B1_D0_a10kst02*B1_kst2*B1_Vm1;
2306
2307 }
2308 else if(g2[i]==1122){//1121
2309 if(B1_D0_a10kst0==-1){
2310 B1_D0_a10kst01 = BarrierN(1,sD0,sa102,skst1,rDs2,mD0M);
2311 B1_D0_a10kst02 = BarrierN(1,sD0,sa101,skst2,rDs2,mD0M);
2312 }
2313
2314 if(B1_a10_kst==-1){
2315 B1_kst1=BarrierN(1,skst1,sKs,sPi01,rRes2,mass2[i]);
2316 B1_kst2=BarrierN(1,skst2,sKs,sPi02,rRes2,mass2[i]);
2317 }
2318 //B1_Vp1 B1_Vp2--rho+
2319 if(B2_a10_rho==-1){
2320 B1_a11m=BarrierN(2,sa101,srhom1,sPip,rRes2,mass1[i]);
2321
2322 B2_a12m=BarrierN(2,sa102,srhom2,sPip,rRes2,mass1[i]);
2323 }
2324
2325 calt1(pa101,pkst2,t1D0_a101);
2326 calt1(pa102,pkst1,t1D0_a102);
2327 //t1rhop1 t1rhop2
2328 calt2(prhom1,Pip,t1_a102);
2329 calt2(prhom2,Pip,t2_a102);
2330 for(int i=0; i<4; i++)
2331 {
2332 double t1_a10_rhom1 = t1rhom1[i]*G[i][i];
2333 double t1_a10_rhom2 = t1rhom2[i]*G[i][i];
2334
2335 for(int j=0; j<4; j++)
2336 {
2337 double t1_a10_kst01 = t1kst1[j]*G[j][j];
2338 double t1_a10_kst02 = t1kst2[j]*G[j][j];
2339 for(int k=0; k<4; k++)
2340 {
2341 double t1_D0_a101kst02 = t1D0_a101[k]*G[k][k];
2342 double t1_D0_a102kst01 = t1D0_a102[k]*G[k][k];
2343
2344 for(int l=0; l<4; l++)
2345 {
2346 double pD0_a10kst0 = pD0[l]*G[l][l];
2347
2348 for(int m=0; m<4; m++)
2349 {
2350 temp_PDF1 += t1_a10_rhom2*t1_a10_kst01*t1_D0_a101kst02*pD0_a10kst0*t2_a102[m][j]*G[m][m]*E[i][k][l][m];
2351 temp_PDF2 += t1_a10_rhom1*t1_a10_kst02*t1_D0_a102kst01*pD0_a10kst0*t1_a102[m][j]*G[m][m]*E[i][k][l][m];
2352 }
2353 }
2354 }
2355
2356 }
2357 }
2358 temp_PDF1 = temp_PDF1*B1_D0_a10kst01*B1_kst1*B1_Vm2*B2_a12m;
2359
2360 temp_PDF2 = temp_PDF2*B1_D0_a10kst02*B1_kst2*B1_Vm1*B1_a11m;
2361
2362 }
2363 else if(g2[i]==1220){//2101
2364 if(B2_D0_a10kst0==-1){
2365 B2_D0_a10kst01 = BarrierN(2,sD0,sa102,skst1,rDs2,mD0M);
2366 B2_D0_a10kst02 = BarrierN(2,sD0,sa101,skst2,rDs2,mD0M);
2367 }
2368 if(B1_a10_kst==-1){
2369 B1_kst1=BarrierN(1,skst1,sKs,sPi01,rRes2,mass2[i]);
2370 B1_kst2=BarrierN(1,skst2,sKs,sPi02,rRes2,mass2[i]);
2371 }
2372 //B1_Vp1 B1_Vp2--rho+
2373
2374 //t1rhop1 t1rhop2
2375 calt2(pa101,pkst2,t2D0_a101);
2376 calt2(pa102,pkst1,t2D0_a102);
2377
2378 calG2(prhom1,Pip,p1a1m_1);
2379 calG2(prhom2,Pip,p1a1m_2);
2380 for(int i=0; i<4;i++)
2381 {
2382 double t1_a10_rhom1 = t1rhom1[i]*G[i][i];
2383 double t1_a10_rhom2 = t1rhom2[i]*G[i][i];
2384
2385 for(int j=0; j<4; j++)
2386 {
2387 double t1_a10_kst01 = t1kst1[j]*G[j][j];
2388 double t1_a10_kst02 = t1kst2[j]*G[j][j];
2389
2390 for(int k=0; k<4; k++)
2391 {
2392 temp_PDF1 += t1_a10_rhom2*t1_a10_kst01*t2D0_a102[k][i]*p1a1m_2[k][j];
2393 temp_PDF2 += t1_a10_rhom1*t1_a10_kst02*t2D0_a101[k][i]*p1a1m_1[k][j];
2394 }
2395 }
2396 }
2397
2398 temp_PDF1 = temp_PDF1*B2_D0_a10kst01*B1_kst1*B1_Vm2;
2399 temp_PDF2 = temp_PDF2*B2_D0_a10kst02*B1_kst2*B1_Vm1;
2400
2401 }
2402 else if(g2[i]==1222){//2121
2403 if(B2_D0_a10kst0==-1){
2404 B2_D0_a10kst01 = BarrierN(2,sD0,sa102,skst1,rDs2,mD0M);
2405 B2_D0_a10kst02 = BarrierN(2,sD0,sa101,skst2,rDs2,mD0M);
2406 }
2407 if(B1_a10_kst==-1){
2408 B1_kst1=BarrierN(1,skst1,sKs,sPi01,rRes2,mass2[i]);
2409 B1_kst2=BarrierN(1,skst2,sKs,sPi02,rRes2,mass2[i]);
2410 }
2411 if(B2_a10_rho==-1){
2412 B1_a11m=BarrierN(2,sa101,srhom1,sPip,rRes2,mass1[i]);
2413
2414 B2_a12m=BarrierN(2,sa102,srhom2,sPip,rRes2,mass1[i]);
2415 }
2416 //B1_Vp1 B1_Vp2--rho+
2417
2418
2419 calt2(pa101,pkst2,t2D0_a101);
2420 calt2(pa102,pkst1,t2D0_a102);
2421
2422 calt2(prhom1,Pip,t1_a103);
2423 calt2(prhom2,Pip,t2_a103);
2424 //t1rhop1 t1rhop2
2425
2426 for(int i=0; i<4;i++)
2427 {
2428 double t1_a10_rhom1 = t1rhom1[i]*G[i][i];
2429 double t1_a10_rhom2 = t1rhom2[i]*G[i][i];
2430 for(int j=0; j<4; j++)
2431 {
2432 double t1_a10_kst01 = t1kst1[j]*G[j][j];
2433 double t1_a10_kst02 = t1kst2[j]*G[j][j];
2434
2435 for(int k=0; k<4; k++)
2436 {
2437 temp_PDF1 += t1_a10_rhom2*t1_a10_kst01*t2D0_a102[k][i]*t2_a103[k][j]*G[k][k];
2438 temp_PDF2 += t1_a10_rhom1*t1_a10_kst02*t2D0_a101[k][i]*t1_a103[k][j]*G[k][k];
2439 }
2440 }
2441 }
2442
2443 temp_PDF1 = temp_PDF1*B2_D0_a10kst01*B1_kst1*B1_Vm2*B2_a12m;
2444 temp_PDF2 = temp_PDF2*B2_D0_a10kst02*B1_kst2*B1_Vm1*B1_a11m;
2445
2446 }
2447
2448
2449 if(g0[i] == 1){//a1
2450 if(g2[i]==1020||g2[i]==1120||g2[i]==1220) {
2451 propagatorRBW(mass1sq,mass1[i],width1[i],sa101,srhom1,sPip,rRes2,0,proa_m1);
2452 propagatorRBW(mass1sq,mass1[i],width1[i],sa102,srhom2,sPip,rRes2,0,proa_m2);
2453 // printf("shuchu1: %lf, %lf\n", proa_m1[0],proa_m2[0]);
2454 }
2455 else{
2456 propagatorRBW(mass1sq,mass1[i],width1[i],sa101,srhom1,sPip,rRes2,2,proa_m1);
2457 propagatorRBW(mass1sq,mass1[i],width1[i],sa102,srhom2,sPip,rRes2,2,proa_m2);
2458 }
2459
2460 }
2461 else if(g0[i]==0){
2462 proa_m1[0] = 1; proa_m1[1] = 0;
2463 proa_m2[0] = 1; proa_m2[1] = 0;
2464
2465 }
2466
2467 //kst
2468 if(g1[i]==1){
2469 propagatorRBW(mass2sq,mass2[i],width2[i],skst1,sKs,sPi01,rRes2,1,pro1_1);//kst1
2470 propagatorRBW(mass2sq,mass2[i],width2[i],skst2,sKs,sPi02,rRes2,1,pro1_2);//kst2
2471
2472 // printf("shuchu1: %lf, %lf\n", pro1_1[0],pro1_2[0]);
2473 }
2474 else if(g1[i]==0){
2475 pro1_1[0] = 1; pro1_1[1] = 0;
2476 pro1_2[0] = 1; pro1_2[1] = 0;
2477 }
2478
2479 //rho
2480 /* if(r1[i]==1){
2481 propagatorGS(0.60102807,0.77526,0.1478,srhom1,sPi01,sPim,rRes2,pro1Vm1);
2482 propagatorGS(0.60102807,0.77526,0.1478,srhom2,sPi02,sPim,rRes2,pro1Vm2);
2483 }
2484 else if(r1[i]==0){
2485 pro1Vm1[0] = 1; pro1Vm1[1] = 0;
2486 pro1Vm2[0] = 1; pro1Vm2[1] = 0;
2487 }
2488 */
2489
2490 Com_Multi(proa_m1,pro1Vm1,tpro_m1); Com_Multi(tpro_m1,pro1_2,pro_12);//1p
2491 Com_Multi(proa_m2,pro1Vm2,tpro_m2); Com_Multi(tpro_m2,pro1_1,pro_22);//2p
2492 // printf("shuchu2: %lf, %lf\n", amp_tmp1[1],amp_tmp2[1]);
2493
2494 // amp_tmp1[0] =temp_PDF1*pro_22[0];//b1--omega2
2495 // amp_tmp1[1] =temp_PDF1*pro_22[1];
2496 // amp_tmp2[0] =temp_PDF2*pro_12[0];//b1--omega1
2497 // amp_tmp2[1] =temp_PDF2*pro_12[1];
2498
2499 amp_tmp1[0] =-1*(temp_PDF1*pro_22[0]);//b1--omega2
2500 amp_tmp1[1] =-1*(temp_PDF1*pro_22[1]);
2501 amp_tmp2[0] =-1*(temp_PDF2*pro_12[0]);//b1--omega1
2502 amp_tmp2[1] =-1*(temp_PDF2*pro_12[1]);
2503 // printf("shuchu4: %.20f,%.20f,%.20f, %.20f\n", amp_tmp1[1],amp_tmp2[1],amp_tmp1[0],amp_tmp2[0]);
2504 }
2505
2506 else if(modetype[i]==21){
2507 if(g2[i]==0){ //0101---SS
2508 /* if(B1_a1p_kstm_rhop1==-1){
2509 B1_a1p_kstm_rhop1 = BarrierN(1,srhop1,sPip,sPi01,rRes2,mass2[i]);
2510 B1_a1p_kstm_rhop2 = BarrierN(1,srhop2,sPip,sPi02,rRes2,mass2[i]);
2511 }*/
2512
2513 if(B1_kst==-1){
2514 B1_kst01 = BarrierN(1,skst1,sKs,sPi01,rRes2,mass2[i]);
2515 B1_kst02 = BarrierN(1,skst2,sKs,sPi02,rRes2,mass2[i]);
2516
2517 }
2518
2519 if(B1_a10_f0==-1){
2520 B1_a101=BarrierN(1,sa101,sf02,sPi01,rRes2,mass1[i]);//1or2 sf01-->pip pim
2521 B1_a102=BarrierN(1,sa102,sf02,sPi02,rRes2,mass1[i]);
2522 }
2523 //t(a1)
2524 calt1(Pf02,Pi01,t1a101);
2525 calt1(Pf02,Pi02,t1a102);//Pf02-->pip pim
2526
2527 for(int m=0; m<4; m++)
2528 {
2529 double t1_a101=t1a101[m]*G[m][m];
2530 double t1_a102=t1a102[m]*G[m][m];
2531
2532 for(int j=0; j<4; j++)
2533 {
2534 temp_PDF1 += t1kst1[j]*G[j][j]*t1_a102;
2535 temp_PDF2 += t1kst2[j]*G[j][j]*t1_a101;
2536 }
2537
2538 }
2539
2540 temp_PDF1 = temp_PDF1*B1_kst01*B1_a102;
2541 temp_PDF2 = temp_PDF2*B1_kst02*B1_a101;
2542
2543 }
2544 else if(g2[i]==1){
2545 if(B1_D0_a1kst==-1){
2546 B1_D0_a10kst01 = BarrierN(1,sD0,sa101,skst2,rDs2,mD0M);
2547 B1_D0_a10kst02 = BarrierN(1,sD0,sa102,skst1,rDs2,mD0M);
2548 }
2549 if(B1_kst==-1){
2550 B1_kst01 = BarrierN(1,skst1,sKs,sPi01,rRes2,mass2[i]);
2551 B1_kst02 = BarrierN(1,skst2,sKs,sPi02,rRes2,mass2[i]);
2552 }
2553 if(B1_a10_f0==-1){
2554 B1_a101=BarrierN(1,sa101,sf02,sPi01,rRes2,mass1[i]);//1or2 sf01-->pip pim
2555 B1_a102=BarrierN(1,sa102,sf02,sPi02,rRes2,mass1[i]);
2556 }
2557 //B1_Vp1 B1_Vp2--rho+
2558
2559 //Sn
2560 calt1(Pf02,Pi01,t1a101);
2561 calt1(Pf02,Pi02,t1a102);//Pf02-->pip pim
2562
2563
2564 calt1(pa101,pkst2,t1D0a101);
2565 calt1(pa102,pkst1,t1D0a102);
2566
2567 for(int j=0; j<4; j++)
2568 {
2569 double t1_kst1 = t1kst1[j]*G[j][j];
2570 double t1_kst2 = t1kst2[j]*G[j][j];
2571
2572 for(int k=0; k<4; k++)
2573 {
2574 double t1_D0_a101 = t1D0a101[k]*G[k][k];
2575 double t1_D0_a102 = t1D0a102[k]*G[k][k];
2576
2577 for(int l=0; l<4; l++)
2578 {
2579 double pD0_a1pkstm = pD0[l]*G[l][l];
2580
2581 for(int m=0; m<4; m++)
2582 {
2583 temp_PDF1 += t1_kst1*t1_D0_a102*pD0_a1pkstm*t1a102[m]*G[m][m]*E[j][k][l][m];
2584 temp_PDF2 += t1_kst2*t1_D0_a101*pD0_a1pkstm*t1a101[m]*G[m][m]*E[j][k][l][m];
2585 }
2586 }
2587 }
2588 }
2589
2590 temp_PDF1 = temp_PDF1*B1_D0_a10kst02*B1_kst1*B1_a102;
2591 temp_PDF2 = temp_PDF2*B1_D0_a10kst01*B1_kst2*B1_a101;
2592
2593 }
2594 else if(g2[i]==2){//2101
2595 if(B1_D0_a1kst==-1){
2596 B1_D0_a10kst01 = BarrierN(2,sD0,sa101,skst2,rDs2,mD0M);
2597 B1_D0_a10kst02 = BarrierN(2,sD0,sa102,skst1,rDs2,mD0M);
2598 }
2599 if(B1_kst==-1){
2600 B1_kst01 = BarrierN(1,skst1,sKs,sPi01,rRes2,mass2[i]);
2601 B1_kst02 = BarrierN(1,skst2,sKs,sPi02,rRes2,mass2[i]);
2602 }
2603 if(B1_a10_f0==-1){
2604 B1_a101=BarrierN(1,sa101,sf02,sPi01,rRes2,mass1[i]);//1or2 sf01-->pip pim
2605 B1_a102=BarrierN(1,sa102,sf02,sPi02,rRes2,mass1[i]);
2606 }
2607 //B1_Vp1 B1_Vp2--rho+
2608
2609 //Sn
2610 calt1(Pf02,Pi01,t1a101);
2611 calt1(Pf02,Pi02,t1a102);//Pf02-->pip pim
2612
2613
2614 calt2(pa101,pkst2,t2D0a101);
2615 calt2(pa102,pkst1,t2D0a102);
2616
2617 for(int j=0; j<4; j++)
2618 {
2619 double t1_kst1 = t1kst1[j]*G[j][j];
2620 double t1_kst2 = t1kst2[j]*G[j][j];
2621 for(int m=0; m<4; m++)
2622 {
2623 temp_PDF1 += t1_kst1*t2D0a102[m][j]*t1a102[m]*G[m][m];
2624 temp_PDF2 += t1_kst2*t2D0a101[m][j]*t1a101[m]*G[m][m];
2625 }
2626 }
2627
2628 temp_PDF1 = temp_PDF1*B1_D0_a10kst02*B1_kst1*B1_a102;
2629 temp_PDF2 = temp_PDF2*B1_D0_a10kst01*B1_kst2*B1_a101;
2630
2631 }
2632
2633
2634 if(g0[i] == 1){//a1
2635 propagatorRBW(mass1sq,mass1[i],width1[i],sa101,sf02,sPi01,rRes2,1,pro00_01);//sf01-->pi+ pi-
2636 propagatorRBW(mass1sq,mass1[i],width1[i],sa102,sf02,sPi02,rRes2,1,pro00_02);
2637 }
2638 else if(g0[i]==0){
2639 pro00_01[0] = 1; pro00_01[1] = 0;
2640 pro00_02[0] = 1; pro00_02[1] = 0;
2641 }
2642
2643 if(r1[i] == 1){//a1
2644 propagatorsigma500(sf02, sPip, sPim, pf0_2);
2645 }
2646 else if(r1[i]==0){
2647 pf0_2[0] = 1; pf0_2[1] = 0;
2648 }
2649 //kst
2650 if(g1[i]==1){
2651 propagatorRBW(mass2sq,mass2[i],width2[i],skst1,sKs,sPi01,rRes2,1,pro1_3);
2652 propagatorRBW(mass2sq,mass2[i],width2[i],skst2,sKs,sPi02,rRes2,1,pro1_4);
2653 }
2654 else if(g1[i]==0){
2655 pro1_3[0] = 1; pro1_3[1] = 0;
2656 pro1_4[0] = 1; pro1_4[1] = 0;
2657 }
2658
2659
2660 Com_Multi(pro00_01,pro1_4,pfo_a101); Com_Multi(pfo_a101,pf0_2,pro_a101);
2661 Com_Multi(pro00_02,pro1_3,pfo_a102); Com_Multi(pfo_a102,pf0_2,pro_a102);
2662
2663 amp_tmp1[0] = temp_PDF1*pro_a102[0];
2664 amp_tmp1[1] = temp_PDF1*pro_a102[1];
2665 amp_tmp2[0] = temp_PDF2*pro_a101[0];
2666 amp_tmp2[1] = temp_PDF2*pro_a101[1];
2667 }
2668 else if(modetype[i]==34){
2669 if(g2[i]==0){
2670 /* if(B1_D0_a1pkpiSW==-1){
2671 B1_D0_a1pkpiSW = BarrierN(1,sD0,sa1p,skstm2,rDs2,mD0M);
2672 }*/
2673 if(B1_kstm==-1){
2674 B1_kstm = BarrierN(1,skst1,sKs,sPim,rRes2,mass2[i]);
2675
2676 }
2677 if(B1_a10_f0==-1){
2678 B1_a10_f0 = BarrierN(1,sa1p,sf01,sPip,rRes2,mass1[i]);
2679 }
2680 //t(a1)
2681 calt1(Pf01,Pip,t1a1p);//pi01 pi02
2682
2683 // calt1(pa1p,pkstm2,t1D0_a1pkpimSW);
2684 calt1(PKs,Pim,t1kstm);
2685 for(int m=0; m<4; m++)
2686 {
2687 double t1_a1p=t1a1p[m]*G[m][m];
2688
2689 for(int j=0; j<4; j++)
2690 {
2691 temp_PDF1 += t1kstm[j]*G[j][j]*t1_a1p;
2692 }
2693
2694 }
2695 temp_PDF1 = temp_PDF1*B1_kstm*B1_a10_f0;
2696
2697 }
2698 /* else if(g2[i]==1){
2699 if(B1_D0_a1kst==-1){
2700 B1_D0_a1pkstm = BarrierN(1,sD0,sa1p,skstm,rDs2,mD0M);
2701 }
2702 if(B1_kst==-1){
2703 B1_kstm = BarrierN(1,skstm,sKs,sPim,rRes2,mass2[i]); //m2==kst Fn(kst)
2704 }
2705 if(B1_a10_f0==-1){
2706 B1_a1p=BarrierN(1,sa1p,sf01,sPip,rRes2,mass1[i]);
2707 }
2708
2709 calt1(Pf01,Pip,t1a1p);
2710
2711 calt1(PKs,Pim,t1kstm);
2712
2713 calt1(pa1p,pkstm,t1D0a1p);
2714
2715 for(int j=0; j<4; j++)
2716 {
2717 double t1_kstm = t1kstm[j]*G[j][j];
2718
2719 for(int k=0; k<4; k++)
2720 {
2721 double t1_D0_a1p = t1D0a1p[k]*G[k][k];
2722
2723 for(int l=0; l<4; l++)
2724 {
2725 double pD0_a1pkstm = pD0[l]*G[l][l];
2726
2727 for(int m=0; m<4; m++)
2728 {
2729 temp_PDF1 += t1_kstm*t1_D0_a1p*pD0_a1pkstm*t1a1p[m]*G[m][m]*E[j][k][l][m];
2730 }
2731 }
2732 }
2733 }
2734 temp_PDF1 = temp_PDF1*B1_D0_a1pkstm*B1_kstm*B1_a1p;
2735
2736 }
2737 else if(g2[i]==2){//2101
2738 if(B1_D0_a1kst==-1){
2739 B1_D0_a1pkstm = BarrierN(2,sD0,sa1p,skstm,rDs2,mD0M);
2740 }
2741 if(B1_kst==-1){
2742 B1_kstm = BarrierN(1,skstm,sKs,sPim,rRes2,mass2[i]); //m2==kst Fn(kst)
2743 }
2744 if(B1_a10_f0==-1){
2745 B1_a1p=BarrierN(1,sa1p,sf01,sPip,rRes2,mass1[i]);
2746 }
2747
2748
2749 calt1(Pf01,Pip,t1a1p);
2750
2751 calt1(PKs,Pim,t1kstm);
2752
2753 calt2(pa1p,pkstm,t2D0a1p);
2754
2755 for(int j=0; j<4; j++)
2756 {
2757 double t1_kstm = t1kstm[j]*G[j][j];
2758 for(int m=0; m<4; m++)
2759 {
2760 temp_PDF1 += t1_kstm*t2D0a1p[m][j]*t1a1p[m]*G[m][m];
2761 }
2762 }
2763
2764 temp_PDF1 = temp_PDF1*B1_D0_a1pkstm*B1_kstm*B1_a1p;
2765
2766 }*/
2767
2768
2769 if(g0[i] == 1){//a1
2770 propagatorRBW(mass1sq,mass1[i],width1[i],sa1p,sf01,sPip,rRes2,1,pro0_p1);
2771 }
2772 else if(g0[i]==0){
2773 pro0_p1[0] = 1; pro0_p1[1] = 0;
2774 }
2775 if(r1[i] == 1){//a1
2776 propagatorsigma500(sf01, sPi01, sPi02, pf0_1);
2777 }
2778 else if(r1[i]==0){
2779 pf0_1[0] = 1; pf0_1[1] = 0;
2780 }
2781 //kst
2782 if(g1[i]==1){
2783 // KPiSLASS(skstm,sKs,sPim,pro1_1);
2784 propagatorRBW(mass2sq,mass2[i],width2[i],skstm,sKs,sPim,rRes2,1,pro1_1);
2785 }
2786 else if(g1[i]==0){
2787 pro1_1[0] = 1; pro1_1[1] = 0;
2788 }
2789
2790
2791 Com_Multi(pro0_p1,pro1_1,pfo_a1p); Com_Multi(pfo_a1p,pf0_1,pro_a1p);
2792
2793 amp_tmp1[0] = temp_PDF1*pro_a1p[0];
2794 amp_tmp1[1] = temp_PDF1*pro_a1p[1];
2795 amp_tmp2[0] = 0;
2796 amp_tmp2[1] = 0;
2797 }
2798 else if (modetype[i]==23){
2799 if(g2[i]==0){ //0101---SS
2800 if(B1_D0_f1==-1){
2801 B1_D0_f1 = BarrierN(1,sD0,sf1285,sKs,rDs2,mD0M); //m2==kst Fn(kst)
2802 }
2803
2804 calG2(prho0,Pf01,p2f1285);
2805 calt1(pf1285,PKs,t1D0_f1285);
2806
2807 //t1rhop1 t1rhop2
2808
2809 for(int i=0; i<4; i++)
2810 {
2811 double t1_f1_rho0 = t1rho0[i]*G[i][i];
2812
2813 for(int j=0; j<4; j++)
2814 {
2815 temp_PDF1 += t1_f1_rho0*t1D0_f1285[j]*G[j][j]*p2f1285[i][j];
2816 }
2817
2818 }
2819
2820 temp_PDF1 = temp_PDF1*B1_D0_f1*B1_V0;
2821
2822 }
2823 if (g2[i]==1)//spin=0 D to rho omega
2824 {
2825 if(B1_D0_f1==-1){
2826 B1_D0_f1 = BarrierN(1,sD0,sf1285,sKs,rDs2,mD0M); //m2==kst Fn(kst)
2827 }
2828 if(B1_f1285==-1) {
2829 B1_f1285 = BarrierN(1,sf1285,srho0,sf01,rRes2,mass1[i]);
2830 }
2831
2832 calt1(prho0,Pf01,t1f1285);
2833 calG2(prho0,Pf01,p2f1285);
2834 calt1(pf1285,PKs,t1D0_f1285);
2835
2836 for(int i=0; i<4; i++)
2837 {
2838 for(int j=0; j<4; j++)
2839 {
2840 double tempf1285 = pf1285[j]*G[j][j]*G[i][i];// p(omega)
2841
2842 for(int k=0; k<4; k++)
2843 {
2844 double temt1f1285 = t1f1285[k]*G[k][k];//Sn(ome)
2845
2846 for(int l=0; l<4; l++)
2847 {
2848 if((i==j)||(i==k)||(i==l)||(j==k)||(j==l)||(k==l)) continue;
2849
2850 double temt1rho0 = t1rho0[l]*G[l][l]*E[i][j][k][l];
2851
2852 for(int m=0; m<4; m++)
2853 {
2854 temp_PDF1 += G[m][m]*t1D0_f1285[m]*p2f1285[m][i]*tempf1285*temt1f1285*temt1rho0;
2855 }
2856 }
2857 }
2858 }
2859 }
2860
2861 temp_PDF1 = temp_PDF1*B1_D0_f1*B1_f1285*B1_V0;
2862 }
2863
2864 else if(g2[i]==2){
2865 if(B1_D0_f1==-1){
2866 B1_D0_f1 = BarrierN(1,sD0,sf1285,sKs,rDs2,mD0M); //m2==kst Fn(kst)
2867 }
2868 //P(a1)
2869 calt2(prho0,Pf01,t2f1285);
2870
2871 calt1(pf1285,PKs,t1D0);
2872
2873 //t1rhop1 t1rhop2
2874
2875 for(int i=0; i<4; i++)
2876 {
2877 double t1_f1_rho0 = t1rho0[i]*G[i][i];
2878
2879 for(int j=0; j<4; j++)
2880 {
2881 temp_PDF1 += t1_f1_rho0*t1D0[j]*G[j][j]*t2f1285[i][j];
2882 }
2883
2884 }
2885
2886 temp_PDF1 = temp_PDF1*B1_D0_f1*B1_V0;
2887
2888 }
2889 if(g0[i] == 1){//a1
2890 propagatorRBW(mass1sq,mass1[i],width1[i],sf1285,srho0,sf01,rRes2,g2[i],pro0_p1);
2891 }
2892 else if(g0[i]==0){
2893 pro0_p1[0] = 1; pro0_p1[1] = 0;
2894 }
2895
2896 //rho
2897 if(r1[i]==1){
2898 propagatorGS(mass2sq,mass2[i],width2[i],srho0,sPip,sPim,rRes2,pro1V0);
2899
2900 }
2901 else if(r1[i]==0){
2902 pro1Vp1[0] = 1; pro1Vp1[1] = 0;
2903 }
2904 if(g1[i]==1){
2905 propagatorsigma500(sf01,sPi01,sPi02,pro1_1);
2906 }
2907 else if(g1[i]==0){
2908 pro1_1[0] = 1; pro1_1[1] = 0;
2909 }
2910
2911 Com_Multi(pro0_p1,pro1V0,tpro_p1); Com_Multi(tpro_p1,pro1_1,pro_11);
2912
2913 amp_tmp1[0] = temp_PDF1*pro_11[0];
2914 amp_tmp1[1] = temp_PDF1*pro_11[1];
2915 amp_tmp2[0] = 0;
2916 amp_tmp2[1] = 0;
2917
2918 }
2919 else if (modetype[i]==25){
2920 if(g2[i]==0)
2921 {
2922 if(B1_omega10==-1) {
2923
2924 B1_omega10=BarrierN(1,somega1,srho01,sPi01,rRes2,mass2[i]);
2925 B1_omega1p=BarrierN(1,somega1,srhop1,sPim,rRes2,mass2[i]);
2926 B1_omega1m=BarrierN(1,somega1,srhom1,sPip,rRes2,mass2[i]);
2927
2928 B1_omega20=BarrierN(1,somega2,srho02,sPi02,rRes2,mass2[i]);
2929 B1_omega2p=BarrierN(1,somega2,srhop2,sPim,rRes2,mass2[i]);
2930 B1_omega2m=BarrierN(1,somega2,srhom2,sPip,rRes2,mass2[i]);
2931
2932 }
2933
2934 if(B1_D0_k1400==-1) {
2935
2936 B1_D0_k14001=BarrierN(1,sD0,sk14001,sPi02,rDs2,mD0M);
2937 B1_D0_k14002=BarrierN(1,sD0,sk14002,sPi01,rDs2,mD0M);
2938 }
2939
2940 for(int i=0; i<4; i++)
2941 {
2942 for(int j=0; j<4; j++)
2943 {
2944 double tempomega1 = pomega1[j]*G[j][j]*G[i][i];
2945 double tempomega2 = pomega2[j]*G[j][j]*G[i][i];
2946
2947 for(int k=0; k<4; k++)
2948 {
2949 double temt1omega01 = t1omega01[k]*G[k][k];
2950 double temt1omega02 = t1omega02[k]*G[k][k];
2951 double temt1omegap1 = t1omegap1[k]*G[k][k];
2952 double temt1omegap2 = t1omegap2[k]*G[k][k];
2953 double temt1omegam1 = t1omegam1[k]*G[k][k];
2954 double temt1omegam2 = t1omegam2[k]*G[k][k];
2955
2956 for(int l=0; l<4; l++)
2957 {
2958 double temt1rho0 = t1rho0[l]*G[l][l]*E[i][j][k][l];
2959 double temt1rhom2 = t1rhom2[l]*G[l][l]*E[i][j][k][l];
2960 double temt1rhop1 = t1rhop1[l]*G[l][l]*E[i][j][k][l];
2961 double temt1rhop2 = t1rhop2[l]*G[l][l]*E[i][j][k][l];
2962 double temt1rhom1 = t1rhom1[l]*G[l][l]*E[i][j][k][l];
2963
2964 for(int m=0; m<4; m++)
2965 {
2966 double temt1D0_k14001 = t1D0_k14001[m]*G[m][m];
2967 double temt1D0_k14002 = t1D0_k14002[m]*G[m][m];
2968 for(int n=0; n<4; n++)
2969 {
2970
2971 if((i==j)||(i==k)||(i==l)||(j==k)||(j==l)||(k==l)) continue;
2972 temp_PDF11 += temt1D0_k14002*pk1400_2[m][n]*p1omega02[n][i]*tempomega2*temt1omega02*temt1rho0*G[n][n];
2973 temp_PDF12 += temt1D0_k14002*pk1400_2[m][n]*p1omegap2[n][i]*tempomega2*temt1omegap2*temt1rhop2*G[n][n];
2974 temp_PDF13 += temt1D0_k14002*pk1400_2[m][n]*p1omegam2[n][i]*tempomega2*temt1omegam2*temt1rhom2*G[n][n];
2975 temp_PDF21 += temt1D0_k14001*pk1400_1[m][n]*p1omega01[n][i]*tempomega1*temt1omega01*temt1rho0*G[n][n];
2976 temp_PDF22 += temt1D0_k14001*pk1400_1[m][n]*p1omegap1[n][i]*tempomega1*temt1omegap1*temt1rhop1*G[n][n];
2977 temp_PDF23 += temt1D0_k14001*pk1400_1[m][n]*p1omegam1[n][i]*tempomega1*temt1omegam1*temt1rhom1*G[n][n];
2978 }
2979 }
2980 }
2981 }
2982 }
2983 }
2984 temp_PDF11 = temp_PDF11*B1_D0_k14002*B1_omega20*B1_V02;
2985 temp_PDF12 = temp_PDF12*B1_D0_k14002*B1_omega2p*B1_Vp2;
2986 temp_PDF13 = temp_PDF13*B1_D0_k14002*B1_omega2m*B1_Vm2;
2987 temp_PDF21 = temp_PDF21*B1_D0_k14001*B1_omega10*B1_V01;
2988 temp_PDF22 = temp_PDF22*B1_D0_k14001*B1_omega1p*B1_Vp1;
2989 temp_PDF23 = temp_PDF23*B1_D0_k14001*B1_omega1m*B1_Vm1;
2990 }
2991 else if(g2[i]==2)
2992 {
2993 if(B1_omega10==-1) {
2994
2995 B1_omega10=BarrierN(1,somega1,srho01,sPi01,rRes2,mass2[i]);
2996 B1_omega1p=BarrierN(1,somega1,srhop1,sPim,rRes2,mass2[i]);
2997 B1_omega1m=BarrierN(1,somega1,srhom1,sPip,rRes2,mass2[i]);
2998
2999 B1_omega20=BarrierN(1,somega2,srho02,sPi02,rRes2,mass2[i]);
3000 B1_omega2p=BarrierN(1,somega2,srhop2,sPim,rRes2,mass2[i]);
3001 B1_omega2m=BarrierN(1,somega2,srhom2,sPip,rRes2,mass2[i]);
3002 }
3003
3004 if(B1_D0_k1400==-1) {
3005
3006 B1_D0_k14001=BarrierN(1,sD0,sk14001,sPi02,rDs2,mD0M);
3007 B1_D0_k14002=BarrierN(1,sD0,sk14002,sPi01,rDs2,mD0M);
3008 }
3009 if(B2_k1400==-1) {
3010 B2_k1400_1=BarrierN(2,sk14001,somega1,sKs,rRes2,mass1[i]);
3011 B2_k1400_2=BarrierN(2,sk14002,somega2,sKs,rRes2,mass1[i]);
3012 }
3013 for(int i=0; i<4; i++)
3014 {
3015 for(int j=0; j<4; j++)
3016 {
3017 double tempomega1 = pomega1[j]*G[j][j]*G[i][i];
3018 double tempomega2 = pomega2[j]*G[j][j]*G[i][i];
3019
3020 for(int k=0; k<4; k++)
3021 {
3022 double temt1omega01 = t1omega01[k]*G[k][k];
3023 double temt1omega02 = t1omega02[k]*G[k][k];
3024 double temt1omegap1 = t1omegap1[k]*G[k][k];
3025 double temt1omegap2 = t1omegap2[k]*G[k][k];
3026 double temt1omegam1 = t1omegam1[k]*G[k][k];
3027 double temt1omegam2 = t1omegam2[k]*G[k][k];
3028
3029 for(int l=0; l<4; l++)
3030 {
3031 double temt1rho0 = t1rho0[l]*G[l][l]*E[i][j][k][l];
3032 double temt1rhop1 = t1rhop1[l]*G[l][l]*E[i][j][k][l];
3033 double temt1rhop2 = t1rhop2[l]*G[l][l]*E[i][j][k][l];
3034 double temt1rhom1 = t1rhom1[l]*G[l][l]*E[i][j][k][l];
3035 double temt1rhom2 = t1rhom2[l]*G[l][l]*E[i][j][k][l];
3036
3037 for(int m=0; m<4; m++)
3038 {
3039 double temt1D0_k14001 = t1D0_k14001[m]*G[m][m];
3040 double temt1D0_k14002 = t1D0_k14002[m]*G[m][m];
3041
3042 for(int n=0; n<4; n++)
3043 {
3044 if((i==j)||(i==k)||(i==l)||(j==k)||(j==l)||(k==l)) continue;
3045 temp_PDF11 += temt1D0_k14002*t2k1400_2[m][n]*p1omega02[n][i]*tempomega2*temt1omega02*temt1rho0*G[n][n];
3046 temp_PDF12 += temt1D0_k14002*t2k1400_2[m][n]*p1omegap2[n][i]*tempomega2*temt1omegap2*temt1rhop2*G[n][n];
3047 temp_PDF13 += temt1D0_k14002*t2k1400_2[m][n]*p1omegam2[n][i]*tempomega2*temt1omegam2*temt1rhom2*G[n][n];
3048 temp_PDF21 += temt1D0_k14001*t2k1400_1[m][n]*p1omega01[n][i]*tempomega1*temt1omega01*temt1rho0*G[n][n];
3049 temp_PDF22 += temt1D0_k14001*t2k1400_1[m][n]*p1omegap1[n][i]*tempomega1*temt1omegap1*temt1rhop1*G[n][n];
3050 temp_PDF23 += temt1D0_k14001*t2k1400_1[m][n]*p1omegam1[n][i]*tempomega1*temt1omegam1*temt1rhom1*G[n][n];
3051 }
3052 }
3053 }
3054 }
3055 }
3056 }
3057 temp_PDF11 = temp_PDF11*B1_D0_k14002*B2_k1400_2*B1_omega20*B1_V02;
3058 temp_PDF12 = temp_PDF12*B1_D0_k14002*B2_k1400_2*B1_omega2p*B1_Vp2;
3059 temp_PDF13 = temp_PDF13*B1_D0_k14002*B2_k1400_2*B1_omega2m*B1_Vm2;
3060
3061 temp_PDF21 = temp_PDF21*B1_D0_k14001*B2_k1400_1*B1_omega10*B1_V01;
3062 temp_PDF22 = temp_PDF22*B1_D0_k14001*B2_k1400_1*B1_omega1p*B1_Vp1;
3063 temp_PDF23 = temp_PDF23*B1_D0_k14001*B2_k1400_1*B1_omega1m*B1_Vm1;
3064 }
3065 if(g0[i]==1) { //b1 omega
3066
3067 propagatorRBW(mass1sq,mass1[i],width1[i],sk14001,somega1,sKs,rRes2,g2[i],prob1_1);
3068 propagatorRBW(mass1sq,mass1[i],width1[i],sk14002,somega2,sKs,rRes2,g2[i],prob1_2);
3069
3070 } else if(g0[i]==0) {
3071 prob1_1[0] = 1; prob1_1[1] = 0;
3072 prob1_2[0] = 1; prob1_2[1] = 0;
3073 }
3074 if(g1[i]==1) { // omega--rho
3075
3076 propagatorRBW(mass2sq,mass2[i],width2[i],somega1,srho01,sPi01,rRes2,1,pro1_01);//01
3077 propagatorRBW(mass2sq,mass2[i],width2[i],somega1,srhop1,sPim,rRes2,1,pro1_p1);//p1
3078 propagatorRBW(mass2sq,mass2[i],width2[i],somega1,srhom1,sPip,rRes2,1,pro1_m1);//m1
3079 propagatorRBW(mass2sq,mass2[i],width2[i],somega2,srho02,sPi02,rRes2,1,pro1_02);//02
3080 propagatorRBW(mass2sq,mass2[i],width2[i],somega2,srhop2,sPim,rRes2,1,pro1_p2);//p2
3081 propagatorRBW(mass2sq,mass2[i],width2[i],somega2,srhom2,sPip,rRes2,1,pro1_m2);//m2
3082 } else if(g1[i]==0) {
3083
3084 pro1_01[0] = 1; pro1_01[1] = 0; pro1_p1[0] = 1; pro1_p1[1] = 0; pro1_m1[0] = 1; pro1_m1[1] = 0;
3085 pro1_02[0] = 1; pro1_02[1] = 0; pro1_p2[0] = 1; pro1_p2[1] = 0; pro1_m2[0] = 1; pro1_m2[1] = 0;
3086
3087 }
3088
3089 //rho--pipi
3090 /* if(r1[i]==1){
3091 propagatorGS(0.60102807,0.77526,0.1478,srho01,sPip,sPim,rRes2,pro1V01);
3092 propagatorGS(0.60102807,0.77526,0.1478,srho02,sPip,sPim,rRes2,pro1V02);
3093 propagatorGS(0.60102807,0.77526,0.1478,srhop1,sPip,sPi01,rRes2,pro1Vp1);
3094 propagatorGS(0.60102807,0.77526,0.1478,srhop2,sPip,sPi02,rRes2,pro1Vp2);
3095 propagatorGS(0.60102807,0.77526,0.1478,srhom1,sPi01,sPim,rRes2,pro1Vm1);
3096 propagatorGS(0.60102807,0.77526,0.1478,srhom2,sPi02,sPim,rRes2,pro1Vm2);
3097
3098 }
3099 else if(r1[i]==0){
3100 pro1V01[0] = 1; pro1V01[1] = 0;
3101 pro1V02[0] = 1; pro1V02[1] = 0;
3102 pro1Vp1[0] = 1; pro1Vp1[1] = 0;
3103 pro1Vp2[0] = 1; pro1Vp2[1] = 0;
3104 pro1Vm1[0] = 1; pro1Vm1[1] = 0;
3105 pro1Vm2[0] = 1; pro1Vm2[1] = 0;
3106 }
3107 */
3108 Com_Multi(prob1_1,pro1_01,tpro_01); Com_Multi(tpro_01,pro1V01,pro_21);//10
3109 Com_Multi(prob1_1,pro1_p1,tpro_p1); Com_Multi(tpro_p1,pro1Vp1,pro_22);//1p
3110 Com_Multi(prob1_1,pro1_m1,tpro_m1); Com_Multi(tpro_m1,pro1Vm1,pro_23);//1m
3111
3112 Com_Multi(prob1_2,pro1_02,tpro_02); Com_Multi(tpro_02,pro1V02,pro_11);//20
3113 Com_Multi(prob1_2,pro1_p2,tpro_p2); Com_Multi(tpro_p2,pro1Vp2,pro_12);//2p
3114 Com_Multi(prob1_2,pro1_m2,tpro_m2); Com_Multi(tpro_m2,pro1Vm2,pro_13);//2m
3115
3116 amp_tmp1[0] = temp_PDF11*pro_11[0] + temp_PDF12*pro_12[0] + temp_PDF13*pro_13[0];//b1--omega2
3117 amp_tmp1[1] = temp_PDF11*pro_11[1] + temp_PDF12*pro_12[1] + temp_PDF13*pro_13[1];
3118 amp_tmp2[0] = temp_PDF21*pro_21[0] + temp_PDF22*pro_22[0] + temp_PDF23*pro_23[0];//b1--omega1
3119 amp_tmp2[1] = temp_PDF21*pro_21[1] + temp_PDF22*pro_22[1] + temp_PDF23*pro_23[1];
3120 }
3121 else if(modetype[i]==44){//K1270 to ks rho+
3122 if(g2[i]==0){ //0101---SS
3123
3124 calG2(prhop1,PKs,p1_k127p1);
3125 calG2(prhop2,PKs,p1_k127p2);
3126 calt1(Pim,Pi01,t1rhom111);
3127 calt1(Pim,Pi02,t1rhom112);
3128 for(int i=0; i<4; i++)
3129 {
3130 double temt1rhom1 = t1rhom111[i]*G[i][i];
3131 double temt1rhom2 = t1rhom112[i]*G[i][i];
3132 for(int j=0; j<4; j++)
3133 {
3134 temp_PDF1 += temt1rhom2*t1rhop1[j]*G[j][j]*p1_k127p1[i][j];
3135 temp_PDF2 += temt1rhom1*t1rhop2[j]*G[j][j]*p1_k127p2[i][j];
3136 }
3137 }
3138
3139 temp_PDF1 = temp_PDF1*B1_Vp1*B1_Vm112;
3140 temp_PDF2 = temp_PDF2*B1_Vp2*B1_Vm111;
3141 }
3142 if(g0[i] == 1){//a1
3143 propagatorRBW(mass1sq,mass1[i],width1[i],sk127p1,srhop1,sKs,rRes2,0,proa_p1);
3144 propagatorRBW(mass1sq,mass1[i],width1[i],sk127p2,srhop2,sKs,rRes2,0,proa_p2);
3145 }
3146 else if(g0[i]==0){
3147 proa_p1[0] = 1; proa_p1[1] = 0;
3148 proa_p2[0] = 1; proa_p2[1] = 0;
3149 }
3150
3151
3152 if(g1[i]==1){
3153 propagatorGS(mass2sq,mass2[i],width2[i],srhom111,sPim,sPi01,rRes2,pro111Vm1);
3154 propagatorGS(mass2sq,mass2[i],width2[i],srhom112,sPim,sPi02,rRes2,pro111Vm2);
3155 }
3156 else if(g1[i]==0){
3157 pro111Vm1[0] = 1; pro111Vm1[1] = 0;
3158 pro111Vm2[0] = 1; pro111Vm2[1] = 0;
3159 }
3160
3161
3162 Com_Multi(proa_p2,pro1Vp2,tpro_p1); Com_Multi(tpro_p1,pro111Vm1,pro_12);//1p
3163 Com_Multi(proa_p1,pro1Vp1,tpro_p2); Com_Multi(tpro_p2,pro111Vm2,pro_22);//2p
3164
3165 amp_tmp1[0] = temp_PDF1*pro_22[0] ;//b1--omega2
3166 amp_tmp1[1] = temp_PDF1*pro_22[1] ;
3167 amp_tmp2[0] = temp_PDF2*pro_12[0] ;//b1--omega1
3168 amp_tmp2[1] = temp_PDF2*pro_12[1] ;
3169 }
3170 else if(modetype[i]==45){//K1270 to ks rho+
3171 if(g2[i]==0){ //0101---SS
3172
3173 calG2(prhom11,PKs,p1_k127m1);
3174 calG2(prhom12,PKs,p1_k127m2);
3175 calt1(Pim,Pi01,t1rhom11);
3176 calt1(Pim,Pi02,t1rhom12);
3177 for(int i=0; i<4; i++)
3178 {
3179 double temt1rhop1 = t1rhop1[i]*G[i][i];
3180 double temt1rhop2 = t1rhop2[i]*G[i][i];
3181 for(int j=0; j<4; j++)
3182 {
3183 temp_PDF1 += temt1rhop2*t1rhom11[j]*G[j][j]*p1_k127m1[i][j];
3184 temp_PDF2 += temt1rhop1*t1rhom12[j]*G[j][j]*p1_k127m2[i][j];
3185 }
3186 }
3187
3188 temp_PDF1 = temp_PDF1*B1_Vp2*B1_Vm11;
3189 temp_PDF2 = temp_PDF2*B1_Vp1*B1_Vm12;
3190 }
3191 if(g0[i] == 1){//a1
3192 propagatorRBW(mass1sq,mass1[i],width1[i],sk127m1,srhom11,sKs,rRes2,0,proa_p1);
3193 propagatorRBW(mass1sq,mass1[i],width1[i],sk127m2,srhom12,sKs,rRes2,0,proa_p2);
3194 }
3195 else if(g0[i]==0){
3196 proa_p1[0] = 1; proa_p1[1] = 0;
3197 proa_p2[0] = 1; proa_p2[1] = 0;
3198 }
3199
3200
3201 if(g1[i]==1){
3202 propagatorGS(mass2sq,mass2[i],width2[i],srhop1,sPip,sPi01,rRes2,pro11Vp1);
3203 propagatorGS(mass2sq,mass2[i],width2[i],srhop2,sPip,sPi02,rRes2,pro11Vp2);
3204 }
3205 else if(g1[i]==0){
3206 pro1_1[0] = 1; pro1_1[1] = 0;
3207 pro1_2[0] = 1; pro1_2[1] = 0;
3208 }
3209
3210
3211 Com_Multi(proa_p2,pro11Vm2,tpro_p1); Com_Multi(tpro_p1,pro11Vp1,pro_12);//1p
3212 Com_Multi(proa_p1,pro11Vm1,tpro_p2); Com_Multi(tpro_p2,pro11Vp2,pro_22);//2p
3213
3214 amp_tmp1[0] = temp_PDF1*pro_22[0] ;//b1--omega2
3215 amp_tmp1[1] = temp_PDF1*pro_22[1] ;
3216 amp_tmp2[0] = temp_PDF2*pro_12[0] ;//b1--omega1
3217 amp_tmp2[1] = temp_PDF2*pro_12[1] ;
3218 }
3219 else if(modetype[i]==46){
3220 if(g2[i]==1){
3221
3222 if(B1_D0_k1270==-1){
3223 B1_D0_k1270 = BarrierN(1,sD0,sk1270,sf01,rDs2,mD0M); //m2==kst Fn(kst)
3224 }
3225 calG2(prho0,PKs,p2a1270);
3226 calt1(pk1270,Pf01,t1D0);
3227 for(int i=0; i<4; i++)
3228 {
3229 double t1_a1270_rho0 = t1rho0[i]*G[i][i];
3230
3231 for(int j=0; j<4; j++)
3232 {
3233 temp_PDF1 += t1_a1270_rho0*t1D0[j]*G[j][j]*p2a1270[i][j];
3234 }
3235
3236 }
3237
3238 temp_PDF1 = temp_PDF1*B1_D0_k1270*B1_V0;
3239 }
3240 else if(g2[i]==2){
3241
3242 if(B1_D0_k1270==-1){
3243 B1_D0_k1270 = BarrierN(1,sD0,sk1270,sf01,rDs2,mD0M); //m2==kst Fn(kst)
3244 }
3245
3246 if(B2_k1270==-1){
3247 B2_k1270 = BarrierN(2,sk1270,sKs,srho0,rRes2,mass1[i]);
3248 }
3249
3250 calt2(PKs,prho0,t2k1270);
3251
3252 calt1(pk1270,Pf01,t1D0);
3253
3254 for(int i=0; i<4; i++)
3255 {
3256 double t1_a1270_rho0 = t1rho0[i]*G[i][i];
3257
3258 for(int j=0; j<4; j++)
3259 {
3260 temp_PDF1 += t1_a1270_rho0*t1D0[j]*G[j][j]*t2k1270[j][i];
3261 }
3262
3263 }
3264 temp_PDF1 = temp_PDF1*B1_D0_k1270*B1_V0*B2_k1270;
3265
3266 }
3267 if(g0[i] == 1){//a1
3268 propagatorRBW(mass1sq,mass1[i],width1[i],sk1270,srho0,sKs,rRes2,g2[i],pro0_p1);
3269 }
3270 else if(g0[i]==0){
3271 pro0_p1[0] = 1; pro0_p1[1] = 0;
3272 }
3273
3274 if(g1[i]==1){
3275 propagatorsigma500(sf01, sPi01, sPi02, pf0_1);
3276 }
3277 else if(g1[i]==0){
3278 pf0_1[0] = 1; pf0_1[1] = 0;
3279 }
3280
3281 if(r1[i]==1){
3282 propagatorGS(0.60102807,0.77526,0.1478,srho0,sPip,sPim,rRes2,pro1V0);
3283
3284 }
3285 else if(r1[i]==0){
3286 pro1V0[0] = 1; pro1V0[1] = 0;
3287 }
3288
3289 Com_Multi(pro0_p1,pro1V0,tpro_p1); Com_Multi(tpro_p1,pf0_1,pro_11);
3290
3291 amp_tmp1[0] = temp_PDF1*pro_11[0];
3292 amp_tmp1[1] = temp_PDF1*pro_11[1];
3293 amp_tmp2[0] = 0;
3294 amp_tmp2[1] = 0;
3295 }
3296
3297
3298 else if(modetype[i]==31){
3299 if(g2[i]==1020){ //0101---SS
3300
3301 if(B1_kst==-1){
3302 B1_kst1=BarrierN(1,skst1,sKs,sPi01,rRes2,0.89166);
3303 B1_kst2=BarrierN(1,skst2,sKs,sPi02,rRes2,0.89166);
3304 }
3305 //P(a1)
3306 calG2(pkst1,Pi02,p2k1400_1);
3307 calG2(pkst2,Pi01,p2k1400_2);
3308
3309
3310 //t1rhop1 t1rhop2
3311
3312 for(int i=0; i<4; i++)
3313 {
3314 double t1_k1400_rho0 = t1rho0[i]*G[i][i];
3315
3316 for(int j=0; j<4; j++)
3317 {
3318 temp_PDF1 += t1_k1400_rho0*t1kst1[j]*G[j][j]*p2k1400_1[i][j];
3319 temp_PDF2 += t1_k1400_rho0*t1kst2[j]*G[j][j]*p2k1400_2[i][j];
3320 }
3321
3322 }
3323
3324 temp_PDF1 = temp_PDF1*B1_kst1*B1_V0;
3325 temp_PDF2 = temp_PDF2*B1_kst2*B1_V0;
3326
3327 }
3328 else if(g2[i]==1120){
3329 if(B1_D0_k1400rho0==-1){
3330 B1_D0_k1400rho0 = BarrierN(1,sD0,sk1400,srho0,rDs2,mD0M);
3331 }
3332
3333 if(B1_kst==-1){
3334 B1_kst1=BarrierN(1,skst1,sKs,sPi01,rRes2,0.89166);
3335 B1_kst2=BarrierN(1,skst2,sKs,sPi02,rRes2,0.89166);
3336 }
3337
3338 //B1_Vp1 B1_Vp2--rho+
3339
3340 //Sn
3341 calt1(pk1400,prho0,t1D0_k1400);
3342
3343 //P
3344 calG2(pkst1,Pi02,p2k1400_1);
3345 calG2(pkst2,Pi01,p2k1400_2);
3346 for(int i=0; i<4; i++)
3347 {
3348 double t1_k1400_rho0 = t1rho0[i]*G[i][i];
3349
3350 for(int j=0; j<4; j++)
3351 {
3352 double t1_k1400_kst1 = t1kst1[j]*G[j][j];
3353 double t1_k1400_kst2 = t1kst2[j]*G[j][j];
3354
3355 for(int k=0; k<4; k++)
3356 {
3357 double t1_D0_k1400 = t1D0_k1400[k]*G[k][k];
3358
3359 for(int l=0; l<4; l++)
3360 { if((i==j)||(i==k)||(i==l)||(j==k)||(j==l)||(k==l)) continue;
3361 double pD0_k1400rho0 = pD0[l]*G[l][l];
3362
3363 for(int m=0; m<4; m++)
3364 {
3365 temp_PDF1 += t1_k1400_rho0*t1_k1400_kst1*t1_D0_k1400*pD0_k1400rho0*p2k1400_1[m][j]*G[m][m]*E[i][k][l][m];
3366 temp_PDF2 += t1_k1400_rho0*t1_k1400_kst2*t1_D0_k1400*pD0_k1400rho0*p2k1400_2[m][j]*G[m][m]*E[i][k][l][m];
3367 }
3368 }
3369 }
3370 }
3371 }
3372 temp_PDF1 = temp_PDF1*B1_D0_k1400rho0*B1_kst1*B1_V0;
3373 temp_PDF2 = temp_PDF2*B1_D0_k1400rho0*B1_kst2*B1_V0;
3374
3375 }
3376 else if(g2[i]==1220){//2101
3377 if(B2_D0_k1400rho0==-1){
3378 B2_D0_k1400rho0 = BarrierN(2,sD0,sk1400,srho0,rDs2,mD0M);
3379 }
3380
3381 if(B1_kst==-1){
3382 B1_kst1=BarrierN(1,skst1,sKs,sPi01,rRes2,0.89166);
3383 B1_kst2=BarrierN(1,skst2,sKs,sPi02,rRes2,0.89166);
3384 }
3385
3386 //B1_Vp1 B1_Vp2--rho+
3387
3388 calt2(pk1400,prho0,t2D0_k1400);
3389 calG2(pkst1,Pi02,p2k1400_1);
3390 calG2(pkst2,Pi01,p2k1400_2);
3391
3392 for(int i=0; i<4;i++)
3393 {
3394 double t1_k1400_rho0= t1rho0[i]*G[i][i];
3395
3396 for(int j=0; j<4; j++)
3397 {
3398 double t1_k1400_kst1 = t1kst1[j]*G[j][j];
3399 double t1_k1400_kst2 = t1kst2[j]*G[j][j];
3400 for(int k=0; k<4; k++)
3401 {
3402 temp_PDF1 += t1_k1400_rho0*t1_k1400_kst1*t2D0_k1400[k][i]*p2k1400_1[k][j];
3403 temp_PDF2 += t1_k1400_rho0*t1_k1400_kst2*t2D0_k1400[k][i]*p2k1400_2[k][j];
3404 }
3405 }
3406 }
3407
3408 temp_PDF1 = temp_PDF1*B2_D0_k1400rho0*B1_kst1*B1_V0;
3409 temp_PDF2 = temp_PDF2*B2_D0_k1400rho0*B1_kst2*B1_V0;
3410
3411 }
3412 if(g0[i] == 1){//a1
3413 propagatorRBW(mass1sq,mass1[i],width1[i],sk1400,skst1,sPi02,rRes2,0,pro0_p1);
3414 propagatorRBW(mass1sq,mass1[i],width1[i],sk1400,skst2,sPi01,rRes2,0,pro0_p2);
3415 }
3416 else if(g0[i]==0){
3417 pro0_p1[0] = 1; pro0_p1[1] = 0;
3418 pro0_p2[0] = 1; pro0_p2[1] = 0;
3419 }
3420
3421
3422 //kst
3423 if(g1[i]==1){
3424 propagatorRBW(0.79505756,0.89166,0.0508,skst1,sKs,sPi01,rRes2,1,pro1Vp1);
3425 propagatorRBW(0.79505756,0.89166,0.0508,skst2,sKs,sPi02,rRes2,1,pro1Vp2);
3426 }
3427 else if(g1[i]==0){
3428 pro1Vp1[0] = 1; pro1Vp1[1] = 0;
3429 pro1Vp2[0] = 1; pro1Vp2[1] = 0;
3430 }
3431
3432 //rho
3433 if(r1[i]==1){
3434 propagatorGS(mass2sq,mass2[i],width2[i],srho0,sPip,sPim,rRes2,pro1_1);
3435
3436 }
3437 else if(r1[i]==0){
3438 pro1_1[0] = 1; pro1_1[1] = 0;
3439 }
3440
3441
3442 Com_Multi(pro0_p1,pro1Vp1,tpro_p1); Com_Multi(tpro_p1,pro1_1,pro_11);
3443 Com_Multi(pro0_p2,pro1Vp2,tpro_p2); Com_Multi(tpro_p2,pro1_1,pro_12);
3444
3445 amp_tmp1[0] = temp_PDF1*pro_11[0];
3446 amp_tmp1[1] = temp_PDF1*pro_11[1];
3447 amp_tmp2[0] = temp_PDF2*pro_12[0];
3448 amp_tmp2[1] = temp_PDF2*pro_12[1];
3449
3450 }
3451
3452 else if(modetype[i]==39){
3453 if (g2[i]==1)
3454 {
3455 if(B1_omega10==-1) {
3456
3457 B1_omega10=BarrierN(1,somega1,srho01,sPi01,rRes2,mass2[i]);
3458 B1_omega1p=BarrierN(1,somega1,srhop1,sPim,rRes2,mass2[i]);
3459 B1_omega1m=BarrierN(1,somega1,srhom1,sPip,rRes2,mass2[i]);
3460
3461 B1_omega20=BarrierN(1,somega2,srho02,sPi02,rRes2,mass2[i]);
3462 B1_omega2p=BarrierN(1,somega2,srhop2,sPim,rRes2,mass2[i]);
3463 B1_omega2m=BarrierN(1,somega2,srhom2,sPip,rRes2,mass2[i]);
3464 }
3465
3466 if(B1_rho14501==-1) {
3467
3468 B1_rho14501=BarrierN(1,srho1450,somega1,sPi02,rRes2,mass1[i]);
3469 B1_rho14502=BarrierN(1,srho1450,somega2,sPi01,rRes2,mass1[i]);
3470 }
3471
3472 if(B1_D0_omega1==-1) {
3473
3474 B1_D0_omega1=BarrierN(1,sD0,srho1450,sKs,rDs2,mD0M);
3475 }
3476
3477 calt1(prho1450,PKs,t1D0_rho1450ks0);
3478 calG2(pomega1,Pi02,p1rho14501);calG2(pomega2,Pi01,p1rho14502);
3479 calt1(pomega1,Pi02,t1rho14501);calt1(pomega2,Pi01,t1rho14502);
3480
3481 for(int i=0; i<4; i++)
3482 {
3483 double temt1D0_rho1450ks0 = t1D0_rho1450ks0[i]*G[i][i];
3484
3485 for(int j=0; j<4; j++)
3486 {
3487 double temp1rho14501 = p1rho14501[i][j]*G[j][j];
3488 double temp1rho14502 = p1rho14502[i][j]*G[j][j];
3489
3490
3491 for(int k=0; k<4; k++)
3492 {
3493
3494 for(int l=0; l<4; l++)
3495 {
3496 double temt1t1rho14501 = t1rho14501[l]*G[l][l];
3497 double temt1t1rho14502 = t1rho14502[l]*G[l][l];
3498
3499 for(int m=0; m<4; m++)
3500 { double temprho1450 = prho1450[m]*G[m][m]*E[j][k][l][m];;
3501
3502 for(int n=0; n<4; n++)
3503 {
3504 double tempomega1 = pomega1[n]*G[n][n];
3505 double tempomega2 = pomega2[n]*G[n][n];
3506
3507 for(int o=0; o<4; o++)
3508 {
3509
3510 double temp1omega01 = p1omega01[k][o]*G[k][k]*G[o][o];
3511 double temp1omega02 = p1omega02[k][o]*G[k][k]*G[o][o];
3512 double temp1omegap1 = p1omegap1[k][o]*G[k][k]*G[o][o];
3513 double temp1omegap2 = p1omegap2[k][o]*G[k][k]*G[o][o];
3514 double temp1omegam1 = p1omegam1[k][o]*G[k][k]*G[o][o];
3515 double temp1omegam2 = p1omegam2[k][o]*G[k][k]*G[o][o];
3516 for(int p=0; p<4; p++)
3517 {
3518
3519 double temt1omega01 = t1omega01[p]*G[p][p];
3520 double temt1omega02 = t1omega02[p]*G[p][p];
3521 double temt1omegap1 = t1omegap1[p]*G[p][p];
3522 double temt1omegap2 = t1omegap2[p]*G[p][p];
3523 double temt1omegam1 = t1omegam1[p]*G[p][p];
3524 double temt1omegam2 = t1omegam2[p]*G[p][p];
3525 for(int q=0; q<4; q++)
3526 {
3527 if((j==k)||(j==l)||(j==m)||(k==l)||(k==m)||(l==m)) continue;
3528 if((n==o)||(n==p)||(n==q)||(o==p)||(o==q)||(p==q)) continue;
3529 temp_PDF11 += temprho1450*temt1D0_rho1450ks0*temt1t1rho14502*temp1rho14502*temp1omega02*E[n][o][p][q]*tempomega2*temt1omega02*t1rho01[q]*G[q][q];
3530 temp_PDF12 += temprho1450*temt1D0_rho1450ks0*temt1t1rho14502*temp1rho14502*temp1omegap2*E[n][o][p][q]*tempomega2*temt1omegap2*t1rhop2[q]*G[q][q];
3531 temp_PDF13 += temprho1450*temt1D0_rho1450ks0*temt1t1rho14502*temp1rho14502*temp1omegam2*E[n][o][p][q]*tempomega2*temt1omegam2*t1rhom2[q]*G[q][q];
3532 temp_PDF21 += temprho1450*temt1D0_rho1450ks0*temt1t1rho14501*temp1rho14501*temp1omega01*E[n][o][p][q]*tempomega1*temt1omega01*t1rho02[q]*G[q][q];
3533 temp_PDF22 += temprho1450*temt1D0_rho1450ks0*temt1t1rho14501*temp1rho14501*temp1omegap1*E[n][o][p][q]*tempomega1*temt1omegap1*t1rhop1[q]*G[q][q];
3534 temp_PDF23 += temprho1450*temt1D0_rho1450ks0*temt1t1rho14501*temp1rho14501*temp1omegam1*E[n][o][p][q]*tempomega1*temt1omegam1*t1rhom1[q]*G[q][q];
3535
3536 }
3537 }
3538 }
3539 }
3540 }
3541 }
3542 }
3543 }
3544 }
3545 temp_PDF11 = temp_PDF11*B1_D0_omega2*B1_omega20*B1_rho14502*B1_V02;
3546 temp_PDF12 = temp_PDF12*B1_D0_omega2*B1_omega2p*B1_rho14502*B1_Vp2;
3547 temp_PDF13 = temp_PDF13*B1_D0_omega2*B1_omega2m*B1_rho14502*B1_Vm2;
3548 temp_PDF21 = temp_PDF21*B1_D0_omega1*B1_omega10*B1_rho14501*B1_V01;
3549 temp_PDF22 = temp_PDF22*B1_D0_omega1*B1_omega1p*B1_rho14501*B1_Vp1;
3550 temp_PDF23 = temp_PDF23*B1_D0_omega1*B1_omega1m*B1_rho14501*B1_Vm1;
3551 }
3552 if(g0[i]==1){
3553 propagatorGS(mass1sq,mass1[i],width1[i],srho1450,somega1,sPi02,rRes2,pro1_1);//kst1
3554 propagatorGS(mass1sq,mass1[i],width1[i],srho1450,somega2,sPi01,rRes2,pro1_2);//kst2
3555
3556 }
3557 else if(g0[i]==0){
3558 pro1_1[0] = 1; pro1_1[1] = 0;
3559 pro1_2[0] = 1; pro1_2[1] = 0;
3560 }
3561 //kst
3562 if(g1[i]==1){
3563 propagatorRBW(mass2sq,mass2[i],width2[i],somega1,srho01,sPi01,rRes2,1,pro0_01);//ome1--rho0 pi01
3564 propagatorRBW(mass2sq,mass2[i],width2[i],somega1,srhop1,sPim,rRes2,1,pro0_p1);//ome1--rhop (pi01)
3565 propagatorRBW(mass2sq,mass2[i],width2[i],somega1,srhom1,sPip,rRes2,1,pro0_m1);//ome1--rhom (pi01)
3566 propagatorRBW(mass2sq,mass2[i],width2[i],somega2,srho02,sPi02,rRes2,1,pro0_02);//ome2--rho0 pi02
3567 propagatorRBW(mass2sq,mass2[i],width2[i],somega2,srhop2,sPim,rRes2,1,pro0_p2);//ome2--rhop (pi02)
3568 propagatorRBW(mass2sq,mass2[i],width2[i],somega2,srhom2,sPip,rRes2,1,pro0_m2);//ome2--rhom(pi02)
3569 }
3570 else if(g1[i]==0){
3571 pro0_01[0] = 1; pro0_01[1] = 0; pro0_p1[0] = 1; pro0_p1[1] = 0; pro0_m1[0] = 1; pro0_m1[1] = 0;
3572 pro0_02[0] = 1; pro0_02[1] = 0; pro0_p2[0] = 1; pro0_p2[1] = 0; pro0_m2[0] = 1; pro0_m2[1] = 0;
3573 }
3574 //rho
3575 /* if(r1[i]==1){
3576 propagatorGS(0.60102807,0.77526,0.1478,srho01,sPip,sPim,rRes2,pro1V01);
3577 propagatorGS(0.60102807,0.77526,0.1478,srho02,sPip,sPim,rRes2,pro1V02);
3578 propagatorGS(0.60102807,0.77526,0.1478,srhop1,sPip,sPi01,rRes2,pro1Vp1);
3579 propagatorGS(0.60102807,0.77526,0.1478,srhop2,sPip,sPi02,rRes2,pro1Vp2);
3580 propagatorGS(0.60102807,0.77526,0.1478,srhom1,sPi01,sPim,rRes2,pro1Vm1);
3581 propagatorGS(0.60102807,0.77526,0.1478,srhom2,sPi02,sPim,rRes2,pro1Vm2);
3582
3583 }
3584 else if(r1[i]==0){
3585 pro1V01[0] = 1; pro1V01[1] = 0;
3586 pro1V02[0] = 1; pro1V02[1] = 0;
3587 pro1Vp1[0] = 1; pro1Vp1[1] = 0;
3588 pro1Vp2[0] = 1; pro1Vp2[1] = 0;
3589 pro1Vm1[0] = 1; pro1Vm1[1] = 0;
3590 pro1Vm2[0] = 1; pro1Vm2[1] = 0;
3591 }
3592 */
3593 Com_Multi(pro0_01,pro1V01,tpro_01); Com_Multi(tpro_01,pro1_2,pro_21);
3594 Com_Multi(pro0_p1,pro1Vp1,tpro_p1); Com_Multi(tpro_p1,pro1_2,pro_22);
3595 Com_Multi(pro0_m1,pro1Vm1,tpro_m1); Com_Multi(tpro_m1,pro1_2,pro_23);
3596
3597 Com_Multi(pro0_02,pro1V02,tpro_02); Com_Multi(tpro_02,pro1_1,pro_11);
3598 Com_Multi(pro0_p2,pro1Vp2,tpro_p2); Com_Multi(tpro_p2,pro1_1,pro_12);
3599 Com_Multi(pro0_m2,pro1Vm2,tpro_m2); Com_Multi(tpro_m2,pro1_1,pro_13);
3600
3601 amp_tmp1[0] =-1*(temp_PDF11*pro_11[0]) + temp_PDF12*pro_12[0] + temp_PDF13*pro_13[0];
3602 amp_tmp1[1] =-1*(temp_PDF11*pro_11[1]) + temp_PDF12*pro_12[1] + temp_PDF13*pro_13[1];
3603 amp_tmp2[0] =-1*(temp_PDF21*pro_21[0]) + temp_PDF22*pro_22[0] + temp_PDF23*pro_23[0];
3604 amp_tmp2[1] =-1*(temp_PDF21*pro_21[1]) + temp_PDF22*pro_22[1] + temp_PDF23*pro_23[1];
3605
3606 // amp_tmp1[0] =temp_PDF11*pro_21[0] + temp_PDF12*pro_22[0] + temp_PDF13*pro_23[0];
3607 // amp_tmp1[1] =temp_PDF11*pro_21[1] + temp_PDF12*pro_22[1] + temp_PDF13*pro_23[1];
3608 // amp_tmp2[0] =temp_PDF21*pro_11[0] + temp_PDF22*pro_12[0] + temp_PDF23*pro_13[0];
3609 // amp_tmp2[1] =temp_PDF21*pro_11[1] + temp_PDF22*pro_12[1] + temp_PDF23*pro_13[1];
3610
3611
3612 }
3613
3614 else if(modetype[i]==41){
3615 amp_tmp1[0] = 1.;
3616 amp_tmp1[1] = 0.;
3617 amp_tmp2[0] = 1.;
3618 amp_tmp2[1] = 0.;
3619 }
3620
3621 //////////////////////////////////////////////
3622 amp_tmp[0] = amp_tmp1[0]+amp_tmp2[0];
3623 amp_tmp[1] = amp_tmp1[1]+amp_tmp2[1];
3624 Com_Multi(amp_tmp,cof,amp_PDF);
3625 PDF[0] += amp_PDF[0];
3626 PDF[1] += amp_PDF[1];
3627
3628 }
3629
3630 double value = PDF[0]*PDF[0] + PDF[1]*PDF[1];
3631
3632 if(value <=0) value = 1e-20;
3633 Result = value;
3634}
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")
n
const Int_t n
Definition: DataBase/tau_mode.c:65
g1
TF1 * g1
Definition: DataBase/tau_mode.c:66
exp
EvtComplex exp(const EvtComplex &c)
Definition: EvtComplex.hh:252
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
EvtD0ToKSpipipi0pi0::clone
EvtDecayBase * clone()
Definition: EvtD0ToKSpipipi0pi0.cc:41
EvtD0ToKSpipipi0pi0::getName
void getName(std::string &name)
Definition: EvtD0ToKSpipipi0pi0.cc:37
EvtD0ToKSpipipi0pi0::decay
void decay(EvtParticle *p)
Definition: EvtD0ToKSpipipi0pi0.cc:184
EvtDecayBase::checkSpinDaughter
void checkSpinDaughter(int d1, EvtSpinType::spintype sp)
Definition: EvtDecayBase.cc:533
EvtDecayBase::checkSpinParent
void checkSpinParent(EvtSpinType::spintype sp)
Definition: EvtDecayBase.cc:520
EvtDecayBase::setProbMax
void setProbMax(double prbmx)
Definition: EvtDecayBase.cc:297
EvtDecayBase::checkNDaug
void checkNDaug(int d1, int d2=-1)
Definition: EvtDecayBase.cc:504
EvtDecayBase::getDaugs
EvtId * getDaugs()
Definition: EvtDecayBase.hh:65
EvtDecayBase::checkNArg
void checkNArg(int a1, int a2=-1, int a3=-1, int a4=-1)
Definition: EvtDecayBase.cc:482
EvtDecayProb::setProb
void setProb(double prob)
Definition: EvtDecayProb.hh:34
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
b
const double b
Definition: slope.cxx:9