BOSS 7.1.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDToKmPipPipPi0.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2// Environment:
3// This software is part of models developed at BES collaboration
4// based on the EvtGen framework. If you use all or part
5// of it, please give an appropriate acknowledgement.
6//
7// Copyright Information: See EvtGen/BesCopyright
8// Copyright (A) 2006 Ping Rong-Gang @IHEP
9//
10// Module: EvtDToKmPipPipPi0.cc
11// see https://indico.ihep.ac.cn/event/16632/contributions/49303/attachments/23510/26646/Xin_Zeng_KPiPiPi0.pdf
12//
13// Modification history:
14//
15// Liaoyuan Dong Oct. 01, 2022 Module created
16//------------------------------------------------------------------------
17#include "EvtGenBase/EvtPatches.hh"
18#include <fstream>
19#include <stdlib.h>
20#include "EvtGenBase/EvtParticle.hh"
21#include "EvtGenBase/EvtGenKine.hh"
22#include "EvtGenBase/EvtPDL.hh"
23#include "EvtGenBase/EvtReport.hh"
24#include "EvtGenBase/EvtResonance.hh"
25#include "EvtGenBase/EvtResonance2.hh"
26#include "EvtGenModels/EvtDToKmPipPipPi0.hh"
27#include <string>
28#include "EvtGenBase/EvtConst.hh"
29#include "EvtGenBase/EvtFlatte.hh"
30#include "EvtGenBase/EvtDecayTable.hh"
31#include "TMatrix.h"
32#include "TMatrixD.h"
33#include "TComplex.h"
34
35using std::endl;
36
37EvtDToKmPipPipPi0::~EvtDToKmPipPipPi0() {}
38
39void EvtDToKmPipPipPi0::getName(std::string& model_name){
40 model_name="DToKmPipPipPi0";
41}
42
43EvtDecayBase* EvtDToKmPipPipPi0::clone(){
44 return new EvtDToKmPipPipPi0;
45}
46
47void EvtDToKmPipPipPi0::init(){
48 // check that there are 0 arguments
49 checkNArg(0);
50 checkNDaug(4);
51
52 checkSpinParent(EvtSpinType::SCALAR);
53 checkSpinDaughter(0,EvtSpinType::SCALAR);
54 checkSpinDaughter(1,EvtSpinType::SCALAR);
55 checkSpinDaughter(2,EvtSpinType::SCALAR);
56 checkSpinDaughter(3,EvtSpinType::SCALAR);
57
58//------------------parameters-----------------
59 rho[0] = 1.00; phi[0] = 0.00;
60 rho[1] = 1.7545e-01; phi[1] = 1.4102e+00;
61 rho[2] = -5.1947e-01; phi[2] = 3.1246e+00;
62 rho[3] = 9.9483e-01; phi[3] = 4.6896e-01;
63 rho[4] = 7.1106e-01; phi[4] = 4.3110e-01;
64 rho[5] = 6.3661e+00; phi[5] = 4.9169e+00;
65 rho[6] = 1.4064e+01; phi[6] = -2.3110e+00;
66 rho[7] = 3.0275e+00; phi[7] = 2.9577e+00;
67 rho[8] = -6.9322e+00; phi[8] = 4.6902e+00;
68 rho[9] = 1.1069e+00; phi[9] = 3.4885e-01;
69 rho[10]= 5.5398e+00; phi[10]= -8.4782e-01;
70 rho[11]= -7.5436e-01; phi[11]= 1.0600e+00;
71 rho[12]= 7.1533e-01; phi[12]= -5.3175e+00;
72 rho[13]= 1.7440e+00; phi[13]= -5.1549e+00;
73 mD = 1.86486;
74
75 //cout << "Initializing DToKpipipi0" << endl;
76 //for (int i=0; i<14; i++) {
77 // cout << i << " rho,phi = " << rho[i] << ", "<< phi[i] << endl;
78 //}
79
80//------------------------new----------------------
81 mKst0 = 0.89555;
82 mrho = 0.77511;
83 mrho1450 = 1.465;
84 mK1270 = 1.289;
85 mK1400 = 1.403;
86 mK1460 = 1.482;
87 mK1650 = 1.672;
88 mK1680 = 1.718;
89
90 GKst0 = 0.0473;
91 Grho = 0.1491;
92 Grho1450 = 0.400;
93 GK1270 = 0.116;
94 GK1400 = 0.174;
95 GK1460 = 0.3356;
96 GK1650 = 0.158;
97 GK1680 = 0.322;
98
99 mass_Pion = 0.13957;
100 mass_Pion_N = 0.134977;
101 mass_Eta = 0.547862;
102 mass_Kaon = 0.493677;
103 math_pi = 3.1415926;
104
105 rD2 = 5.0; // 5*5
106 rRes1 = 3.0; // 3*3
107 rRes2 = 3.0; // 3*3
108
109 GS1 = 0.636619783;
110 GS2 = 0.01860182466;
111 GS3 = 0.1591549458; // 1/(2*math_2pi)
112 GS4 = 0.00620060822; // mass_Pion2/math_pi
113
114 int GG[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
115 int EE[4][4][4][4] =
116 { { {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
117 {{0,0,0,0}, {0,0,0,0}, {0,0,0,1}, {0,0,-1,0}},
118 {{0,0,0,0}, {0,0,0,-1}, {0,0,0,0}, {0,1,0,0} },
119 {{0,0,0,0}, {0,0,1,0}, {0,-1,0,0}, {0,0,0,0} } },
120 { {{0,0,0,0}, {0,0,0,0}, {0,0,0,-1}, {0,0,1,0} },
121 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
122 {{0,0,0,1}, {0,0,0,0}, {0,0,0,0}, {-1,0,0,0}},
123 {{0,0,-1,0}, {0,0,0,0}, {1,0,0,0}, {0,0,0,0} } },
124 { {{0,0,0,0}, {0,0,0,1}, {0,0,0,0}, {0,-1,0,0}},
125 {{0,0,0,-1}, {0,0,0,0}, {0,0,0,0}, {1,0,0,0} },
126 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
127 {{0,1,0,0}, {-1,0,0,0}, {0,0,0,0}, {0,0,0,0} } },
128 { {{0,0,0,0}, {0,0,-1,0}, {0,1,0,0}, {0,0,0,0} },
129 {{0,0,1,0}, {0,0,0,0}, {-1,0,0,0}, {0,0,0,0} },
130 {{0,-1,0,0}, {1,0,0,0}, {0,0,0,0}, {0,0,0,0} },
131 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} } } };
132 for (int i=0; i<4; i++) {
133 for (int j=0; j<4; j++) {
134 G[i][j] = GG[i][j];
135 for (int k=0; k<4; k++) {
136 for (int l=0; l<4; l++) {
137 E[i][j][k][l] = EE[i][j][k][l];
138 }
139 }
140 }
141 }
142}
143
144void EvtDToKmPipPipPi0::initProbMax() {
145 setProbMax(100000);
146}
147
148void EvtDToKmPipPipPi0::decay( EvtParticle *p){
149//--------------max
150/*
151 double maxprob=0,maxq1270;
152 for(int ir=0;ir<=20000000;ir++){
153 p->initializePhaseSpace(getNDaug(),getDaugs());
154 EvtVector4R km = p->getDaug(0)->getP4();
155 EvtVector4R pip1 = p->getDaug(1)->getP4();
156 EvtVector4R pip2 = p->getDaug(2)->getP4();
157 EvtVector4R pi0 = p->getDaug(3)->getP4();
158
159 double Pip1[4],Pip2[4],Km[4],Pi0[4];
160 Km[0] = km.get(0); Pip1[0] = pip1.get(0); Pip2[0] = pip2.get(0);Pi0[0] = pi0.get(0);
161 Km[1] = km.get(1); Pip1[1] = pip1.get(1); Pip2[1] = pip2.get(1);Pi0[1] = pi0.get(1);
162 Km[2] = km.get(2); Pip1[2] = pip1.get(2); Pip2[2] = pip2.get(2);Pi0[2] = pi0.get(2);
163 Km[3] = km.get(3); Pip1[3] = pip1.get(3); Pip2[3] = pip2.get(3);Pi0[3] = pi0.get(3);
164
165
166 double _prob;
167 double value;
168 int nstates=14;
169 int modetype[14]= {1,1,8,4,12,13,12,17,8,4,13,9,4,11};
170 int g0[14] = {1,1,1,1,1,0,0,1,1,1,0,1,1,1};
171 int g1[14] = {1,1,1,1,1,1,1,1,0,0,0,0,1,1};
172 int g2[14] = {0,1,0,0,1,1,1,0,2,0,1,0,2,1};
173 double mass1[14] = {mKst0,mKst0,mrho, mKst0, mKst0, mrho, mKst0, mKst0, mrho, mKst0, mrho, mKst0, mKst0, mKst0};
174 double mass2[14] = {mrho, mrho, mK1270,mK1400,mK1460,mK1460,mK1460,mrho, mK1270,mK1400,mK1460,mK1270,mK1400, mK1680};
175 double width1[14] = {GKst0,GKst0,Grho, GKst0, GKst0, Grho, GKst0, GKst0, Grho, GKst0, Grho, GKst0, GKst0, GKst0};
176 double width2[14] = {Grho, Grho, GK1270,GK1400,GK1460,GK1460,GK1460,Grho, GK1270,GK1400,GK1460,GK1270,GK1400, GK1680};
177
178 calEvaMy(Km,Pip1,Pip2,Pi0,mass1,mass2,width1,width2,rho,phi,g0,g1,g2,modetype,nstates,_prob);
179 if(_prob>maxprob) {
180 maxprob=_prob;
181 }
182 }
183 printf("maxprob = %.10f\n", maxprob);
184*/
185//--------------------main--------------
186
187 p->initializePhaseSpace(getNDaug(),getDaugs());
188 EvtVector4R km = p->getDaug(0)->getP4();
189 EvtVector4R pip1 = p->getDaug(1)->getP4();
190 EvtVector4R pip2 = p->getDaug(2)->getP4();
191 EvtVector4R pi0 = p->getDaug(3)->getP4();
192 int num = EvtPDL::getStdHep(p->getId());
193
194 double Pip1[4],Pip2[4],Km[4],Pi0[4];
195 //----------------------------------------------------------------------------------------
196 if(num>0){
197 Km[0] = km.get(0); Pip1[0] = pip1.get(0); Pip2[0] = pip2.get(0);Pi0[0] = pi0.get(0);
198 Km[1] = km.get(1); Pip1[1] = pip1.get(1); Pip2[1] = pip2.get(1);Pi0[1] = pi0.get(1);
199 Km[2] = km.get(2); Pip1[2] = pip1.get(2); Pip2[2] = pip2.get(2);Pi0[2] = pi0.get(2);
200 Km[3] = km.get(3); Pip1[3] = pip1.get(3); Pip2[3] = pip2.get(3);Pi0[3] = pi0.get(3);
201}
202 if(num<0){
203 Km[0] = km.get(0); Pip1[0] = pip1.get(0); Pip2[0] = pip2.get(0);Pi0[0] = pi0.get(0);
204 Km[1] = -km.get(1); Pip1[1] = -pip1.get(1); Pip2[1] = -pip2.get(1);Pi0[1] = -pi0.get(1);
205 Km[2] = -km.get(2); Pip1[2] = -pip1.get(2); Pip2[2] = -pip2.get(2);Pi0[2] = -pi0.get(2);
206 Km[3] = -km.get(3); Pip1[3] = -pip1.get(3); Pip2[3] = -pip2.get(3);Pi0[3] = -pi0.get(3);
207}
208 double value;
209 int nstates=14;
210 int modetype[14]= {1,1,8,4,12,13,12,17,8,4,13,9,4,11};
211 int g0[14] = {1,1,1,1,1,0,0,1,1,1,0,1,1,1};
212 int g1[14] = {1,1,1,1,1,1,1,1,0,0,0,0,1,1};
213 int g2[14] = {0,1,0,0,1,1,1,0,2,0,1,0,2,1};
214 double mass1[14] = {mKst0,mKst0,mrho, mKst0, mKst0, mrho, mKst0, mKst0, mrho, mKst0, mrho, mKst0, mKst0, mKst0};
215 double mass2[14] = {mrho, mrho, mK1270,mK1400,mK1460,mK1460,mK1460,mrho, mK1270,mK1400,mK1460,mK1270,mK1400, mK1680};
216 double width1[14] = {GKst0,GKst0,Grho, GKst0, GKst0, Grho, GKst0, GKst0, Grho, GKst0, Grho, GKst0, GKst0, GKst0};
217 double width2[14] = {Grho, Grho, GK1270,GK1400,GK1460,GK1460,GK1460,Grho, GK1270,GK1400,GK1460,GK1270,GK1400, GK1680};
218
219 calEvaMy(Km,Pip1,Pip2,Pi0,mass1,mass2,width1,width2,rho,phi,g0,g1,g2,modetype,nstates,value);
220 setProb(value);
221
222 return;
223}
224
225void EvtDToKmPipPipPi0::KPiSLASS(double sa, double sb, double sc, double prop[2]) {
226 const double m1430 = 1.441;
227 const double sa0 = 2.076481; // m1430*m1430;
228 const double w1430 = 0.193;
229 const double Lass1 = 0.25/sa0;
230 double tmp = sb-sc;
231 double tmp1 = sa0+tmp;
232 double q0 = Lass1*tmp1*tmp1-sb;
233 if(q0<0) q0 = 1e-16;
234 double tmp2 = sa+tmp;
235 double qs = 0.25*tmp2*tmp2/sa-sb;
236 double q = sqrt(qs);
237 double width = w1430*q*m1430/sqrt(sa*q0);
238 double temp_R = atan(m1430*width/(sa0-sa));
239 if(temp_R<0) temp_R += math_pi;
240 double deltaR = -1.915 + temp_R;
241 double temp_F = atan(0.226*q/(2.0-3.8194*qs)); // 2.0*1.07 = 2.14; 1.8*1.07 = 1.926
242 if(temp_F<0) temp_F += math_pi;
243 double deltaF = 0.002 + temp_F;
244 double deltaS = deltaR + 2.0*deltaF;
245 double t1 = 0.96*sin(deltaF);
246 double t2 = sin(deltaR);
247 double CF[2], CS[2];
248 CF[0] = cos(deltaF);
249 CF[1] = sin(deltaF);
250 CS[0] = cos(deltaS);
251 CS[1] = sin(deltaS);
252 prop[0] = t1*CF[0] + t2*CS[0];
253 prop[1] = t1*CF[1] + t2*CS[1];
254}
255
256
257void EvtDToKmPipPipPi0::Com_Multi(double a1[2], double a2[2], double res[2]){
258 res[0] = a1[0]*a2[0]-a1[1]*a2[1];
259 res[1] = a1[1]*a2[0]+a1[0]*a2[1];
260}
261void EvtDToKmPipPipPi0::Com_Divide(double a1[2], double a2[2], double res[2]){
262 res[0] = (a1[0]*a2[0]+a1[1]*a2[1])/(a2[0]*a2[0]+a2[1]*a2[1]);
263 res[1] = (a1[1]*a2[0]-a1[0]*a2[1])/(a2[0]*a2[0]+a2[1]*a2[1]);
264}
265double EvtDToKmPipPipPi0::SCADot(double a1[4], double a2[4])
266{
267 double _cal = a1[0]*a2[0]-a1[1]*a2[1]-a1[2]*a2[2]-a1[3]*a2[3];
268 return _cal;
269}
270double EvtDToKmPipPipPi0::Barrier(int l, double sa, double sb, double sc, double r,double mass)
271{
272 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
273// if(q < 0) q = 1e-16;
274 if(q < 0) q = -q;
275
276 double z;
277 z = q*r*r;
278 double sa0;
279 sa0 = mass*mass;
280 double q0 = (sa0+sb-sc)*(sa0+sb-sc)/(4*sa0)-sb;
281// if(q0 < 0) q0 = 1e-16;
282 if(q0 < 0) q0 = -q0;
283 double z0 = q0*r*r;
284 double F = 0.0;
285 if(l == 0) F = 1;
286 if(l == 1) F = sqrt((1+z0)/(1+z));
287 if(l == 2) F = sqrt((9+3*z0+z0*z0)/(9+3*z+z*z));
288 // if(l == 1) F = sqrt((2*r)/(1+z));
289 // if(l == 2) F = sqrt((13*r*r)/(9+3*z+z*z));
290 return F;
291}
292void EvtDToKmPipPipPi0::calt1(double daug1[4], double daug0[4], double t1[4]){
293 double p, pq;
294 double pa[4], qa[4];
295 for(int i=0; i<4; i++){
296 pa[i] = daug1[i] + daug0[i];
297 qa[i] = daug1[i] - daug0[i];
298 }
299 p = SCADot(pa,pa);
300 pq = SCADot(pa,qa);
301 for(int i=0; i<4; i++){
302 t1[i] = qa[i] - pq/p*pa[i];
303 }
304}
305void EvtDToKmPipPipPi0::calt2(double daug1[4], double daug0[4], double t2[4][4]){
306 double p, r;
307 double pa[4], t1[4];
308 calt1(daug1,daug0,t1);
309 r = SCADot(t1,t1)/3.0;
310 for(int i=0; i<4; i++) {
311 pa[i] = daug1[i] + daug0[i];
312 }
313 p = SCADot(pa,pa);
314 for(int i=0; i<4; i++) {
315 for(int j=0; j<4; j++) {
316 t2[i][j] = t1[i]*t1[j] - r*(G[i][j]-pa[i]*pa[j]/p);
317 }
318 }
319}
320
321double EvtDToKmPipPipPi0::wid(double mass2, double mass, double sa, double sb, double sc, double r, int l){
322 double widm = 0.;
323 double m = sqrt(sa);
324 double tmp = sb-sc;
325 double tmp1 = sa+tmp;
326 double q = 0.25*tmp1*tmp1/sa-sb;
327// if(q<0) q = 1e-16;
328 if(q<0) q = -q;
329 double tmp2 = mass2+tmp;
330 double q0 = 0.25*tmp2*tmp2/mass2-sb;
331// if(q0<0) q0 = 1e-16;
332 if(q0<0) q0 = -q0;
333 double z = q*r*r;
334 double z0 = q0*r*r;
335 double t = q/q0;
336 if(l == 0) {widm = sqrt(t)*mass/m;}
337 else if(l == 1) {widm = t*sqrt(t)*mass/m*(1+z0)/(1+z);}
338 else if(l == 2) {widm = t*t*sqrt(t)*mass/m*(9+3*z0+z0*z0)/(9+3*z+z*z);}
339 return widm;
340
341}
342double EvtDToKmPipPipPi0::widl1(double mass2, double mass, double sa, double sb, double sc, double r)
343{
344 double widm = 0.;
345 double m = sqrt(sa);
346 double tmp = sb-sc;
347 double tmp1 = sa+tmp;
348 double q = 0.25*tmp1*tmp1/sa-sb;
349// if(q<0) q = 1e-16;
350 if(q<0) q = -q;
351 double tmp2 = mass2+tmp;
352 double q0 = 0.25*tmp2*tmp2/mass2-sb;
353// if(q0<0) q0 = 1e-16;
354 if(q0<0) q0 = -q0;
355 double z = q*r*r;
356 double z0 = q0*r*r;
357 double F = (1+z0)/(1+z);
358 double t = q/q0;
359 widm = t*sqrt(t)*mass/m*F;
360 return widm;
361
362}
363void EvtDToKmPipPipPi0::propagatorNBW(double mass2, double mass, double width, double sa, double sb, double sc, double r, int l, double prop[2]){
364 double a[2], b[2];
365 a[0] = 1;
366 a[1] = 0;
367 b[0] = mass2-sa;
368 b[1] = -mass*width;
369 Com_Divide(a,b,prop);
370}
371
372void EvtDToKmPipPipPi0::propagatorRBW(double mass2, double mass, double width, double sa, double sb, double sc, double r, int l, double prop[2]){
373 double a[2], b[2];
374 a[0] = 1;
375 a[1] = 0;
376 b[0] = mass2-sa;
377 b[1] = -mass*width*wid(mass2,mass,sa,sb,sc,r,l);
378 Com_Divide(a,b,prop);
379}
380void EvtDToKmPipPipPi0::propagatorRBWl1(double mass2, double mass, double width, double sa, double sb, double sc, double r, double prop[2])
381 //firstly code for propagator_omg, now could use for propagatorRBW that l=1, and will be more quick.
382{
383 double a[2], b[2];
384 a[0] = 1;
385 a[1] = 0;
386 b[0] = mass2-sa;
387 b[1] = -mass*width*widl1(mass2,mass,sa,sb,sc,r);
388 Com_Divide(a,b,prop);
389}
390//------------GS---used by rho----------------------------
391void EvtDToKmPipPipPi0::propagatorGS(double mass2, double mass, double width, double sa, double sb, double sc, double r, double prop[2])
392{
393 double a[2], b[2];
394 double tmp = sb-sc;
395 double tmp1 = sa+tmp;
396 double q2 = 0.25*tmp1*tmp1/sa-sb;
397// if(q2<0) q2 = 1e-16;
398 if(q2<0) q2 = -q2;
399
400 double tmp2 = mass2+tmp;
401 double q02 = 0.25*tmp2*tmp2/mass2-sb;
402// if(q02<0) q02 = 1e-16;
403 if(q02<0) q02 = -q02;
404
405 double q = sqrt(q2);
406 double q0 = sqrt(q02);
407 double m = sqrt(sa);
408 double q03 = q0*q02;
409 double tmp3 = log(mass+2*q0)+1.2760418309; // log(mass_2Pion) = 1.2760418309;
410
411 double h = GS1*q/m*(log(m+2*q)+1.2760418309);
412 double h0 = GS1*q0/mass*tmp3;
413 double dh = h0*(0.125/q02-0.5/mass2)+GS3/mass2;
414 double d = GS2/q02*tmp3+GS3*mass/q0-GS4*mass/q03;
415 double f = mass2/q03*(q2*(h-h0)+(mass2-sa)*q02*dh);
416
417 a[0] = 1.0+d*width/mass;
418 a[1] = 0.0;
419 b[0] = mass2-sa+width*f;
420 b[1] = -mass*width*widl1(mass2,mass,sa,sb,sc,r);
421 Com_Divide(a,b,prop);
422}
423void EvtDToKmPipPipPi0::calEvaMy(
424 double* Km, double* Pip1, double* Pip2, double* Pi0,
425 double* mass1, double* mass2,
426 double* width1, double* width2,
427 double* amp, double* phase,
428 int* g0, int* g1, int* g2,
429 int* modetype, int nstates, double & Result)
430
431{
432
433 double pKstr1[4],pKstr2[4],pKstr3[4],prho1[4],prho2[4],pK11[4],pK12[4],pK13[4],pA1[4],pD[4];
434 double Amp_KPiS1[2], Amp_PiPiS1[2], Amp_KPiS2[2], Amp_PiPiS2[2];
435
436 for(int i=0;i!=4;i++){
437 pKstr1[i]=Km[i]+Pip1[i];
438 pKstr2[i]=Km[i]+Pip2[i];
439 pKstr3[i]=Km[i]+Pi0[i];
440 prho1[i]=Pip1[i]+Pi0[i];
441 prho2[i]=Pip2[i]+Pi0[i];
442 pK11[i]=Km[i]+prho1[i];
443 pK12[i]=Km[i]+prho2[i];
444 pK13[i]=Km[i]+Pip1[i]+Pip2[i];
445 pD[i]=Km[i]+prho1[i]+Pip2[i];
446 }
447 double skstr1,skstr2;
448 double skstr3,srho1,srho2,sk11,sk12,sk13,sD;
449 double skaon,spion1,spion2,spi0;
450 skaon = SCADot(Km,Km);
451 spion1 = SCADot(Pip1,Pip1);
452 spion2 = SCADot(Pip2,Pip2);
453 spi0 = SCADot(Pi0,Pi0);
454 skstr1 = SCADot(pKstr1,pKstr1);
455 skstr2 = SCADot(pKstr2,pKstr2);
456 skstr3 = SCADot(pKstr3,pKstr3);
457 srho1 = SCADot(prho1,prho1);
458 srho2 = SCADot(prho2,prho2);
459 sk11 = SCADot(pK11,pK11);
460 sk12 = SCADot(pK12,pK12);
461 sk13 = SCADot(pK13,pK13);
462 sD = SCADot(pD,pD);
463 double t1rho1[4],t1rho2[4],t1kstr3[4];
464//------------------------
465 double t1kstr1[4],t1kstr2[4], t2k21[4][4], t2k22[4][4], t2k31[4][4], t2k32[4][4], t2k11[4][4], t2k12[4][4], t2k13[4][4], t2k14[4][4];
466 double t2rho1[4][4], t2rho2[4][4], t2kstr3[4][4], t2kstr1[4][4], t2kstr2[4][4];
467 double t1K11[4],t1K12[4],t1K13[4],t1K14[4];
468 double t2kk1[4][4],t2kk2[4][4];
469
470//--------------------------
471 calt1(Pip1,Pi0,t1rho1);
472 calt1(Pip2,Pi0,t1rho2);
473 calt1(Km,Pip1,t1kstr1);
474 calt1(Km,Pip2,t1kstr2);
475 calt1(Km,Pi0 ,t1kstr3);
476 calt1(pKstr1,Pi0 ,t1K11);
477 calt1(pKstr2,Pi0 ,t1K12);
478 calt1(pKstr3,Pip1 ,t1K13);
479 calt1(pKstr3,Pip2 ,t1K14);
480 calt2(Pip2,Pi0,t2rho1);
481 calt2(Pip2,Pi0,t2rho2);
482 calt2(Km,Pi0,t2kstr3);
483 calt2(Km,Pip1,t2kstr1);
484 calt2(Km,Pip2,t2kstr2);
485 calt2(pKstr3,Pip1,t2k11);
486 calt2(pKstr3,Pip2,t2k12);
487 calt2(prho1,Km,t2k21);
488 calt2(prho2,Km,t2k22);
489 calt2(pKstr1,Pi0,t2k31);
490 calt2(pKstr2,Pi0,t2k32);
491 calt2(pKstr1,Pip2,t2kk1);
492 calt2(pKstr2,Pip1,t2kk2);
493//----------------------------------
494 double cof[2],amp_tmp1[2],amp_PDF[2], PDF[2];
495 double amp_tmp2[2];
496 double amp_tmp[2];
497
498 PDF[0] = 0.0;
499 PDF[1] = 0.0;
500
501 double temp_PDF1, temp_PDF2, tt1, tt2, tmp1, tmp2;
502
503 double temp_PDF3, temp_PDF4;
504
505 double pro_1[2], pro_2[2], pro0_1[2], pro1_1[2], pro0_2[2], pro1_2[2], temp_1[2], temp_2[2];
506
507 double t1D[4], t1D_1[4], t1D_2[4];
508
509 double t1K1_1[4], t1K1_2[4],t1K1_3[4],t1K1_4[4];
510
511 double B1_rho1=-1.0, B1_rho2=-1.0;
512 double B1_D_K1_Pi2=-1.0, B1_D_K1_Pi1=-1.0, B2_D_K1_Pi2=-1.0, B2_D_K1_Pi1=-1.0;
513 double B1_K1_rho1=-1.0, B1_K1_rho2=-1.0, B2_K1_rho1=-1.0, B2_K1_rho2=-1.0;
514 double B1_K1_Kstr3_Pi1=-1.0, B1_K1_Kstr3_Pi2=-1.0, B2_K1_Kstr3_Pi1=-1.0, B2_K1_Kstr3_Pi2=-1.0;
515 double B1_K1_Kstr1_Pi0=-1.0, B1_K1_Kstr2_Pi0=-1.0;
516 double B1_kstr3=-1.0;
517 double B2_K1p_kstr1=-1.0,B2_K1p_kstr2=-1.0;
518 double B1_K1p_kstr1=-1.0,B1_K1p_kstr2=-1.0;
519 double B1_kstr1_Pi1=-1.0, B1_kstr2_Pi2=-1.0, B1_D_Kstr2_rho1=-1.0, B1_D_Kstr1_rho2=-1.0,B1_D_K1p_Pi0=-1.0;
520 double B2_K1_Kstr1=-1.0, B2_K1_Kstr2=-1.0, B2_D_Kstr1_rho2=-1.0, B2_D_Kstr2_rho1=-1.0;
521 double B2_rho1=-1.0, B2_rho2=-1.0;
522 double B2_kstr1_Pi1=-1.0, B2_kstr2_Pi2=-1.0, B2_kstr3_Pi0=-1.0;;
523 double mass1sq, mass2sq;
524 int isKst=0, isRho=0, isPiPi_S=0, isKPi_S=0;
525 double proRho1[2], proRho2[2], proPiPi1_S[2], proPiPi2_S[2], proKPi_S[2];
526
527 double tmp3, tmp4;
528 double t2D[4][4], t2D_1[4][4], t2D_2[4][4];
529 double pro2_1[2], pro2_2[2], pro3_1[2], pro3_2[2], pro4_1[2], pro4_2[2];
530 double t2f2_1[4][4],t2f2_2[4][4];
531 double t2G2_1[4][4], t2G2_2[4][4], t2D11_1[4][4],t2D11_2[4][4];
532
533//-------------------------------------------------------------------------
534 for(int i=0; i<nstates; i++){
535 cof[0] = amp[i]*cos(phase[i]);
536 cof[1] = amp[i]*sin(phase[i]);
537 mass1sq = mass1[i]*mass1[i];
538 mass2sq = mass2[i]*mass2[i];
539 temp_PDF1 = 0;
540 temp_PDF2 = 0;
541 temp_PDF3 = 0;
542 temp_PDF4 = 0;
543
544 if(modetype[i]==1)
545 {
546 if (B1_kstr1_Pi1<0.0) B1_kstr1_Pi1 = Barrier(1,skstr1,skaon,spion1,rRes2,mass1[i]);
547 if (B1_kstr2_Pi2<0.0) B1_kstr2_Pi2 = Barrier(1,skstr2,skaon,spion2,rRes2,mass1[i]);
548 if (B1_rho1<0.0) B1_rho1 = Barrier(1,srho1,spion1,spi0,rRes2,mass2[i]);
549 if (B1_rho2<0.0) B1_rho2 = Barrier(1,srho2,spion2,spi0,rRes2,mass2[i]);
550 if (g2[i]==0)
551 {
552 for(int w=0; w<4; w++)
553 {
554 temp_PDF1 += G[w][w]*t1kstr1[w]*t1rho2[w];
555 temp_PDF2 += G[w][w]*t1kstr2[w]*t1rho1[w];
556 }
557 tmp1 = B1_kstr1_Pi1*B1_rho2*temp_PDF1;
558 tmp2 = B1_kstr2_Pi2*B1_rho1*temp_PDF2;
559 }
560 else if (g2[i]==1)
561 {
562 calt1(pKstr1,prho2,t1D_1);
563 calt1(pKstr2,prho1,t1D_2);
564 for(int w=0; w<4; w++)
565 {
566 for(int j=0; j<4; j++)
567 {
568 for(int k=0; k<4; k++)
569 {
570 for(int l=0; l<4; l++)
571 {
572 temp_PDF1 += E[w][j][k][l]*pD[w]*t1D_1[j]*t1kstr1[k]*t1rho2[l]*G[w][w]*G[j][j]*G[k][k]*G[l][l];
573 temp_PDF2 += E[w][j][k][l]*pD[w]*t1D_2[j]*t1kstr2[k]*t1rho1[l]*G[w][w]*G[j][j]*G[k][k]*G[l][l];
574 }
575 }
576 }
577 }
578 if (B1_D_Kstr1_rho2<0.0) B1_D_Kstr1_rho2 = Barrier(1,sD,skstr1,srho2,rD2,mD);
579 if (B1_D_Kstr2_rho1<0.0) B1_D_Kstr2_rho1 = Barrier(1,sD,skstr2,srho1,rD2,mD);
580 tmp1 = B1_kstr1_Pi1*B1_rho2*B1_D_Kstr1_rho2*temp_PDF1;
581 tmp2 = B1_kstr2_Pi2*B1_rho1*B1_D_Kstr2_rho1*temp_PDF2;
582 }
583 else if (g2[i]==2)
584 {
585 calt2(pKstr1,prho2,t2D_1);
586 calt2(pKstr2,prho1,t2D_2);
587 for(int w=0; w<4; w++)
588 {
589 for(int j=0; j<4; j++)
590 {
591 temp_PDF1 += t2D_1[w][j]*t1kstr1[w]*t1rho2[j]*G[w][w]*G[j][j];
592 temp_PDF2 += t2D_2[w][j]*t1kstr2[w]*t1rho1[j]*G[w][w]*G[j][j];
593 }
594 }
595 if (B2_D_Kstr1_rho2<0.0) B2_D_Kstr2_rho1 = Barrier(2,sD,skstr1,srho2,rD2,mD);
596 if (B2_D_Kstr2_rho1<0.0) B2_D_Kstr2_rho1 = Barrier(2,sD,skstr2,srho1,rD2,mD);
597 tmp1 = B1_kstr1_Pi1*B1_rho2*B2_D_Kstr1_rho2*temp_PDF1;
598 tmp2 = B1_kstr2_Pi2*B1_rho1*B2_D_Kstr2_rho1*temp_PDF2;
599 }
600 if(g0[i]==1){
601 propagatorRBW(mass1sq,mass1[i],width1[i],skstr1,skaon,spion1,rRes2,1,pro0_1);
602 propagatorRBW(mass1sq,mass1[i],width1[i],skstr2,skaon,spion2,rRes2,1,pro0_2);
603 }
604 else if(g0[i]==0){
605 pro0_1[0] = 1; pro0_1[1] = 0; pro0_2[0] = 1; pro0_2[1] = 0;
606 }
607 if(g1[i]==1){
608 propagatorGS(mass2sq,mass2[i],width2[i],srho2,spion2,spi0,rRes2,pro1_1);
609 propagatorGS(mass2sq,mass2[i],width2[i],srho1,spion1,spi0,rRes2,pro1_2);
610 }
611 else if(g1[i]==0){
612 pro1_1[0] = 1; pro1_1[1] = 0; pro1_2[0] = 1; pro1_2[1] = 0;
613 }
614 Com_Multi(pro0_1,pro1_1,pro_1);
615 Com_Multi(pro0_2,pro1_2,pro_2);
616 amp_tmp1[0] = tmp1*pro_1[0];
617 amp_tmp1[1] = tmp1*pro_1[1];
618 amp_tmp2[0] = tmp2*pro_2[0];
619 amp_tmp2[1] = tmp2*pro_2[1];
620}
621 else if (modetype[i]==4) {
622 if (B1_D_K1_Pi2<0.0) B1_D_K1_Pi2 = Barrier(1,sD,sk11,spion2,rD2,mass2[i]);
623 if (B1_D_K1_Pi1<0.0) B1_D_K1_Pi1 = Barrier(1,sD,sk12,spion1,rD2,mass2[i]);
624 if (B1_kstr3<0.0) B1_kstr3 = Barrier(1,skstr3,skaon,spi0,rRes2,mass1[i]);
625 if (B1_kstr1_Pi1<0.0) B1_kstr1_Pi1= Barrier(1,skstr1,skaon,spion1,rRes2,mass1[i]);
626 if (B1_kstr2_Pi2<0.0) B1_kstr2_Pi2= Barrier(1,skstr2,skaon,spion2,rRes2,mass1[i]);
627 calt1(pK11,Pip2,t1D_1);
628 calt1(pK12,Pip1,t1D_2);
629 if (g2[i]==0) {
630 for(int w=0; w<4; w++) {
631 for(int j=0; j<4; j++) {
632 tt2 = t1kstr3[j]*G[w][w]*G[j][j];
633 temp_PDF1 += t1D_1[w]*(G[w][j]-pK11[w]*pK11[j]/sk11)*t1kstr1[j]*G[w][w]*G[j][j];
634 temp_PDF2 += t1D_2[w]*(G[w][j]-pK12[w]*pK12[j]/sk12)*t1kstr2[j]*G[w][w]*G[j][j];
635 temp_PDF3 += t1D_1[w]*(G[w][j]-pK11[w]*pK11[j]/sk11)*tt2;
636 temp_PDF4 += t1D_2[w]*(G[w][j]-pK12[w]*pK12[j]/sk12)*tt2;
637 }
638 }
639 tmp1 = B1_kstr1_Pi1*B1_D_K1_Pi2*temp_PDF1;
640 tmp2 = B1_kstr2_Pi2*B1_D_K1_Pi1*temp_PDF2;
641 tmp3 = B1_kstr3*B1_D_K1_Pi2*temp_PDF3;
642 tmp4 = B1_kstr3*B1_D_K1_Pi1*temp_PDF4;
643 }
644 else if (g2[i]==2) {
645 for(int w=0; w<4; w++) {
646 for(int j=0; j<4; j++) {
647 tt2 = t1kstr3[j]*G[w][w]*G[j][j];
648 temp_PDF1 += t1D_1[w]*t2k31[w][j]*t1kstr1[j]*G[w][w]*G[j][j];
649 temp_PDF2 += t1D_2[w]*t2k32[w][j]*t1kstr2[j]*G[w][w]*G[j][j];
650 temp_PDF3 += t1D_1[w]*t2k11[w][j]*tt2;
651 temp_PDF4 += t1D_2[w]*t2k12[w][j]*tt2;
652 }
653 }
654 if (B2_K1_Kstr3_Pi1<0.0) B2_K1_Kstr3_Pi1 = Barrier(2,sk11,skstr3,spion1,rRes2,mass2[i]);
655 if (B2_K1_Kstr3_Pi2<0.0) B2_K1_Kstr3_Pi2 = Barrier(2,sk12,skstr3,spion2,rRes2,mass2[i]);
656 if (B2_K1_Kstr1<0.0) B2_K1_Kstr1 = Barrier(2,sk11,skstr1,spi0,rRes2,mass2[i]);
657 if (B2_K1_Kstr2<0.0) B2_K1_Kstr2 = Barrier(2,sk12,skstr2,spi0,rRes2,mass2[i]);
658 tmp1 = B1_kstr1_Pi1*B2_K1_Kstr1*B1_D_K1_Pi2*temp_PDF1;
659 tmp2 = B1_kstr2_Pi2*B2_K1_Kstr2*B1_D_K1_Pi1*temp_PDF2;
660 tmp3 = B1_kstr3*B2_K1_Kstr3_Pi1*B1_D_K1_Pi2*temp_PDF3;
661 tmp4 = B1_kstr3*B2_K1_Kstr3_Pi2*B1_D_K1_Pi1*temp_PDF4;
662 }
663 propagatorRBW(mass1sq,mass1[i],width1[i],skstr3,skaon,spi0,rRes2,1,pro2_1);
664 propagatorRBW(mass1sq,mass1[i],width1[i],skstr3,skaon,spi0,rRes2,1,pro2_2);
665 propagatorRBW(mass2sq,mass2[i],width2[i],sk11,skstr3,spion1,rRes2,g2[i],pro3_1);
666 propagatorRBW(mass2sq,mass2[i],width2[i],sk12,skstr3,spion2,rRes2,g2[i],pro3_2);
667 if (g0[i]==1) {
668 propagatorRBW(mass1sq,mass1[i],width1[i],skstr1,skaon,spion1,rRes2,1,pro0_1);
669 propagatorRBW(mass1sq,mass1[i],width1[i],skstr2,skaon,spion2,rRes2,1,pro0_2);
670 } else if (g0[i]==0) {
671 pro0_1[0] = 1; pro0_1[1] = 0; pro0_2[0] = 1; pro0_2[1] = 0; pro2_1[0] = 1; pro2_1[1] = 0; pro2_2[0] = 1; pro2_2[1] = 0;
672 }
673 if (g1[i]==1) {
674 propagatorRBW(mass2sq,mass2[i],width2[i],sk11,skstr1,spi0,rRes2,g2[i],pro1_1);
675 propagatorRBW(mass2sq,mass2[i],width2[i],sk12,skstr2,spi0,rRes2,g2[i],pro1_2);
676 } else if (g1[i]==0) {
677 pro1_1[0] = 1; pro1_1[1] = 0; pro1_2[0] = 1; pro1_2[1] = 0; pro3_1[0] = 1; pro3_1[1] = 0; pro3_2[0] = 1; pro3_2[1] = 0;
678
679 }
680 Com_Multi(pro0_1,pro1_1,pro_1);
681 Com_Multi(pro2_1,pro3_1,pro4_1);
682 Com_Multi(pro0_2,pro1_2,pro_2);
683 Com_Multi(pro2_2,pro3_2,pro4_2);
684 amp_tmp1[0] = tmp1*pro_1[0]-tmp3*pro4_1[0];
685 amp_tmp1[1] = tmp1*pro_1[1]-tmp3*pro4_1[1];
686 amp_tmp2[0] = tmp2*pro_2[0]-tmp4*pro4_2[0];
687 amp_tmp2[1] = tmp2*pro_2[1]-tmp4*pro4_2[1];
688 }
689 else if (modetype[i]==8) {
690 if (B1_rho1<0.0) B1_rho1 = Barrier(1,srho1,spion1,spi0,rRes2,mass1[i]);
691 if (B1_rho2<0.0) B1_rho2 = Barrier(1,srho2,spion2,spi0,rRes2,mass1[i]);
692 if (B1_D_K1_Pi2<0.0) B1_D_K1_Pi2 = Barrier(1,sD,sk11,spion2,rD2,mD);
693 if (B1_D_K1_Pi1<0.0) B1_D_K1_Pi1 = Barrier(1,sD,sk12,spion1,rD2,mD);
694 calt1(pK11,Pip2,t1D_1);
695 calt1(pK12,Pip1,t1D_2);
696 if(g2[i]==0) {
697 for(int w=0; w<4; w++) {
698 for(int j=0; j<4; j++) {
699 temp_PDF1 += t1D_1[w]*(G[w][j]-pK11[w]*pK11[j]/sk11)*t1rho1[j]*G[w][w]*G[j][j];
700 temp_PDF2 += t1D_2[w]*(G[w][j]-pK12[w]*pK12[j]/sk12)*t1rho2[j]*G[w][w]*G[j][j];
701 }
702 }
703 tmp1 = B1_rho1*B1_D_K1_Pi2*temp_PDF1;
704 tmp2 = B1_rho2*B1_D_K1_Pi1*temp_PDF2;
705 } else if(g2[i]==2) {
706 for(int w=0; w<4; w++) {
707 for(int j=0; j<4; j++) {
708 temp_PDF1 += t1D_1[w]*t2k21[w][j]*t1rho1[j]*G[w][w]*G[j][j];
709 temp_PDF2 += t1D_2[w]*t2k22[w][j]*t1rho2[j]*G[w][w]*G[j][j];
710 }
711 }
712 if (B2_K1_rho1<0.0) B2_K1_rho1 = Barrier(2,sk11,srho1,skaon,rRes2,mass2[i]);
713 if (B2_K1_rho2<0.0) B2_K1_rho2 = Barrier(2,sk12,srho2,skaon,rRes2,mass2[i]);
714 tmp1 = B1_rho1*B2_K1_rho1*B1_D_K1_Pi2*temp_PDF1;
715 tmp2 = B1_rho2*B2_K1_rho2*B1_D_K1_Pi1*temp_PDF2;
716 }
717 if(g0[i]==1) {
718 propagatorGS(mass1sq,mass1[i],width1[i],srho1,spion1,spi0,rRes2,pro0_1);
719 propagatorGS(mass1sq,mass1[i],width1[i],srho2,spion2,spi0,rRes2,pro0_2);
720 } else if(g0[i]==0) {
721 pro0_1[0] = 1; pro0_1[1] = 0; pro0_2[0] = 1; pro0_2[1] = 0;
722 }
723 if(g1[i]==1) {
724 propagatorRBW(mass2sq,mass2[i],width2[i],sk11,srho1,skaon,rRes2,g2[i],pro1_1);
725 propagatorRBW(mass2sq,mass2[i],width2[i],sk12,srho2,skaon,rRes2,g2[i],pro1_2);
726 } else if(g1[i]==0) {
727 pro1_1[0] = 1; pro1_1[1] = 0; pro1_2[0] = 1; pro1_2[1] = 0;
728 }
729
730 Com_Multi(pro0_1,pro1_1,pro_1);
731 Com_Multi(pro0_2,pro1_2,pro_2);
732 amp_tmp1[0] = tmp1*pro_1[0];
733 amp_tmp1[1] = tmp1*pro_1[1];
734 amp_tmp2[0] = tmp2*pro_2[0];
735 amp_tmp2[1] = tmp2*pro_2[1];
736 }
737 else if (modetype[i]==11) {
738 if (B1_kstr3<0.0) B1_kstr3 = Barrier(1,skstr3,skaon,spi0,rRes2,mass1[i]);
739 if (B1_kstr1_Pi1<0.0) B1_kstr1_Pi1 = Barrier(1,skstr1,skaon,spion1,rRes2,mass1[i]);
740 if (B1_kstr2_Pi2<0.0) B1_kstr2_Pi2 = Barrier(1,skstr2,skaon,spion2,rRes2,mass1[i]);
741 if (B1_D_K1_Pi2<0.0) B1_D_K1_Pi2 = Barrier(1,sD,sk11,spion2,rD2,mD);
742 if (B1_D_K1_Pi1<0.0) B1_D_K1_Pi1 = Barrier(1,sD,sk12,spion1,rD2,mD);
743 for(int w=0; w<4; w++) {
744 for(int j=0; j<4; j++) {
745 for(int k=0; k<4; k++){
746 for(int l=0; l<4; l++){
747 temp_PDF1 += E[w][j][k][l]*pK11[w]*G[w][w]*(pKstr1[j]-Pi0[j])*G[j][j]*Pip2[k]*G[k][k]*(Km[l]-Pip1[l])*G[l][l];
748 temp_PDF2 += E[w][j][k][l]*pK12[w]*G[w][w]*(pKstr2[j]-Pi0[j])*G[j][j]*Pip1[k]*G[k][k]*(Km[l]-Pip2[l])*G[l][l];
749 temp_PDF3 += E[w][j][k][l]*pK11[w]*G[w][w]*(pKstr3[j]-Pip1[j])*G[j][j]*Pip2[k]*G[k][k]*(Km[l]-Pi0[l])*G[l][l];
750 temp_PDF4 += E[w][j][k][l]*pK12[w]*G[w][w]*(pKstr3[j]-Pip2[j])*G[j][j]*Pip1[k]*G[k][k]*(Km[l]-Pi0[l])*G[l][l];
751 }
752 }
753 }
754 }
755 if (B1_K1_Kstr1_Pi0<0.0) B1_K1_Kstr1_Pi0 = Barrier(1,sk11,skstr1,spi0,rRes2,mass2[i]);
756 if (B1_K1_Kstr2_Pi0<0.0) B1_K1_Kstr2_Pi0 = Barrier(1,sk12,skstr2,spi0,rRes2,mass2[i]);
757 if (B1_K1_Kstr3_Pi1<0.0) B1_K1_Kstr3_Pi1 = Barrier(1,sk11,skstr3,spion1,rRes2,mass2[i]);
758 if (B1_K1_Kstr3_Pi2<0.0) B1_K1_Kstr3_Pi1 = Barrier(1,sk12,skstr3,spion2,rRes2,mass2[i]);
759 tmp1 = B1_kstr1_Pi1*B1_K1_Kstr1_Pi0*B1_D_K1_Pi2*temp_PDF1;
760 tmp2 = B1_kstr2_Pi2*B1_K1_Kstr2_Pi0*B1_D_K1_Pi1*temp_PDF2;
761 tmp3 = B1_kstr3*B1_K1_Kstr3_Pi1*B1_D_K1_Pi2*temp_PDF3;
762 tmp4 = B1_kstr3*B1_K1_Kstr3_Pi2*B1_D_K1_Pi1*temp_PDF4;
763 propagatorRBW(mass1sq,mass1[i],width1[i],skstr3,skaon,spi0,rRes2,1,pro2_1);
764 propagatorRBW(mass1sq,mass1[i],width1[i],skstr3,skaon,spi0,rRes2,1,pro2_2);
765 propagatorRBW(mass2sq,mass2[i],width2[i],sk11,skstr3,spion1,rRes2,1,pro3_1);
766 propagatorRBW(mass2sq,mass2[i],width2[i],sk12,skstr3,spion2,rRes2,1,pro3_2);
767 if(g0[i]==1) {
768 propagatorRBW(mass1sq,mass1[i],width1[i],skstr1,skaon,spion1,rRes2,1,pro0_1);
769 propagatorRBW(mass1sq,mass1[i],width1[i],skstr2,skaon,spion2,rRes2,1,pro0_2);
770 } else if(g0[i]==0) {
771 pro0_1[0] = 1; pro0_1[1] = 0; pro0_2[0] = 1; pro0_2[1] = 0; pro2_1[0] = 1; pro2_1[1] = 0; pro2_2[0] = 1; pro2_2[1] = 0;
772 }
773 if(g1[i]==1) {
774 propagatorRBW(mass2sq,mass2[i],width2[i],sk11,skstr1,spi0,rRes2,1,pro1_1);
775 propagatorRBW(mass2sq,mass2[i],width2[i],sk12,skstr2,spi0,rRes2,1,pro1_2);
776 } else if(g1[i]==0) {
777 pro1_1[0] = 1; pro1_1[1] = 0; pro1_2[0] = 1; pro1_2[1] = 0; pro3_1[0] = 1; pro3_1[1] = 0; pro3_2[0] = 1; pro3_2[1] = 0;
778 }
779 Com_Multi(pro0_1,pro1_1,pro_1);
780 Com_Multi(pro2_1,pro3_1,pro4_1);
781 Com_Multi(pro0_2,pro1_2,pro_2);
782 Com_Multi(pro2_2,pro3_2,pro4_2);
783 amp_tmp1[0] = tmp1*pro_1[0]-tmp3*pro4_1[0];
784 amp_tmp1[1] = tmp1*pro_1[1]-tmp3*pro4_1[1];
785 amp_tmp2[0] = tmp2*pro_2[0]-tmp4*pro4_2[0];
786 amp_tmp2[1] = tmp2*pro_2[1]-tmp4*pro4_2[1];
787 }
788
789
790 else if (modetype[i]==12) {
791 if (B1_K1_Kstr3_Pi1<0.0) B1_K1_Kstr3_Pi1 = Barrier(1,sk11,skstr3,spion1,rRes2,mass2[i]);
792 if (B1_K1_Kstr3_Pi2<0.0) B1_K1_Kstr3_Pi2 = Barrier(1,sk12,skstr3,spion2,rRes2,mass2[i]);
793 if (B1_K1_Kstr1_Pi0<0.0) B1_K1_Kstr1_Pi0 = Barrier(1,sk11,skstr1,spi0,rRes2,mass2[i]);
794 if (B1_K1_Kstr2_Pi0<0.0) B1_K1_Kstr2_Pi0 = Barrier(1,sk12,skstr2,spi0,rRes2,mass2[i]);
795 if (B1_kstr3<0.0) B1_kstr3 = Barrier(1,skstr3,skaon,spi0,rRes2,mass1[i]);
796 if (B1_kstr1_Pi1<0.0) B1_kstr1_Pi1= Barrier(1,skstr1,skaon,spion1,rRes2,mass1[i]);
797 if (B1_kstr2_Pi2<0.0) B1_kstr2_Pi2= Barrier(1,skstr2,skaon,spion2,rRes2,mass1[i]);
798 for(int w=0; w<4; w++) {
799 temp_PDF1 += G[w][w]*Pi0[w]*t1kstr1[w];
800 temp_PDF2 += G[w][w]*Pi0[w]*t1kstr2[w];
801 temp_PDF3 += G[w][w]*Pip1[w]*t1kstr3[w];
802 temp_PDF4 += G[w][w]*Pip2[w]*t1kstr3[w];
803 }
804 tmp1 = B1_kstr1_Pi1*B1_K1_Kstr1_Pi0*temp_PDF1;
805 tmp2 = B1_kstr2_Pi2*B1_K1_Kstr2_Pi0*temp_PDF2;
806 tmp3 = B1_kstr3*B1_K1_Kstr3_Pi1*temp_PDF3;
807 tmp4 = B1_kstr3*B1_K1_Kstr3_Pi2*temp_PDF4;
808 propagatorRBW(mass1sq,mass1[i],width1[i],skstr3,skaon,spi0,rRes2,1,pro2_1);
809 propagatorRBW(mass1sq,mass1[i],width1[i],skstr3,skaon,spi0,rRes2,1,pro2_2);
810 propagatorRBW(mass2sq,mass2[i],width2[i],sk11,skstr3,spion1,rRes2,g2[i],pro3_1);
811 propagatorRBW(mass2sq,mass2[i],width2[i],sk12,skstr3,spion2,rRes2,g2[i],pro3_2);
812 if (g0[i]==1) {
813 propagatorRBW(mass1sq,mass1[i],width1[i],skstr1,skaon,spion1,rRes2,1,pro0_1);
814 propagatorRBW(mass1sq,mass1[i],width1[i],skstr2,skaon,spion2,rRes2,1,pro0_2);
815 } else if (g0[i]==0) {
816 pro0_1[0] = 1; pro0_1[1] = 0; pro0_2[0] = 1; pro0_2[1] = 0; pro2_1[0] = 1; pro2_1[1] = 0; pro2_2[0] = 1; pro2_2[1] = 0;
817 }
818 if (g1[i]==1) {
819 propagatorRBW(mass2sq,mass2[i],width2[i],sk11,skstr1,spi0,rRes2,g2[i],pro1_1);
820 propagatorRBW(mass2sq,mass2[i],width2[i],sk12,skstr2,spi0,rRes2,g2[i],pro1_2);
821 } else if (g1[i]==0) {
822 pro1_1[0] = 1; pro1_1[1] = 0; pro1_2[0] = 1; pro1_2[1] = 0; pro3_1[0] = 1; pro3_1[1] = 0; pro3_2[0] = 1; pro3_2[1] = 0;
823
824 }
825 Com_Multi(pro0_1,pro1_1,pro_1);
826 Com_Multi(pro2_1,pro3_1,pro4_1);
827 Com_Multi(pro0_2,pro1_2,pro_2);
828 Com_Multi(pro2_2,pro3_2,pro4_2);
829 amp_tmp1[0] = tmp1*pro_1[0]-tmp3*pro4_1[0];
830 amp_tmp1[1] = tmp1*pro_1[1]-tmp3*pro4_1[1];
831 amp_tmp2[0] = tmp2*pro_2[0]-tmp4*pro4_2[0];
832 amp_tmp2[1] = tmp2*pro_2[1]-tmp4*pro4_2[1];
833
834}
835 if(modetype[i]==13)
836 {
837 if(B1_K1_rho1<0.0) B1_K1_rho1 = Barrier(1,sk11,srho1,skaon,rRes2,mass2[i]);
838 if(B1_K1_rho2<0.0) B1_K1_rho2 = Barrier(1,sk12,srho2,skaon,rRes2,mass2[i]);
839 if (B1_rho1<0.0) B1_rho1 = Barrier(1,srho1,spion1,spi0,rRes2,mass1[i]);
840 if (B1_rho2<0.0) B1_rho2 = Barrier(1,srho2,spion2,spi0,rRes2,mass1[i]);
841 for(int w=0; w<4; w++)
842 {
843 temp_PDF1 += G[w][w]*Km[w]*t1rho1[w];
844 temp_PDF2 += G[w][w]*Km[w]*t1rho2[w];
845 }
846 tmp1 = B1_rho1*B1_K1_rho1*temp_PDF1;
847 tmp2 = B1_rho2*B1_K1_rho2*temp_PDF2;
848 if(g1[i]==1){
849 propagatorRBW(mass2sq,mass2[i],width2[i],sk11,skaon,srho1,rRes2,1,pro0_1);
850 propagatorRBW(mass2sq,mass2[i],width2[i],sk12,skaon,srho2,rRes2,1,pro0_2);
851 }
852 else if(g1[i]==0){
853 pro0_1[0] = 1; pro0_1[1] = 0; pro0_2[0] = 1; pro0_2[1] = 0;
854 }
855 if(g0[i]==1){
856 propagatorGS(mass1sq,mass1[i],width1[i],srho1,spion1,spi0,rRes2,pro1_1);
857 propagatorGS(mass1sq,mass1[i],width1[i],srho2,spion2,spi0,rRes2,pro1_2);
858 }
859 else if(g0[i]==0){
860 pro1_1[0] = 1; pro1_1[1] = 0; pro1_2[0] = 1; pro1_2[1] = 0;
861 }
862 Com_Multi(pro0_1,pro1_1,pro_1);
863 Com_Multi(pro0_2,pro1_2,pro_2);
864 amp_tmp1[0] = tmp1*pro_1[0];
865 amp_tmp1[1] = tmp1*pro_1[1];
866 amp_tmp2[0] = tmp2*pro_2[0];
867 amp_tmp2[1] = tmp2*pro_2[1];
868 }
869 else if (modetype[i]==17)
870 {
871 if (B1_rho1<0.0) B1_rho1 = Barrier(1,srho1,spion1,spi0,rRes2,mass2[i]);
872 if (B1_rho2<0.0) B1_rho2 = Barrier(1,srho2,spion2,spi0,rRes2,mass2[i]);
873 if (B1_D_Kstr2_rho1<0.0) B1_D_Kstr2_rho1 = Barrier(1,sD,skstr2,srho1,rD2,mD);
874 if (B1_D_Kstr1_rho2<0.0) B1_D_Kstr1_rho2 = Barrier(1,sD,skstr1,srho2,rD2,mD);
875 calt1(pKstr2,prho1,t1D_1);
876 calt1(pKstr1,prho2,t1D_2);
877 for(int w=0; w<4; w++){
878 temp_PDF1 += t1D_1[w]*t1rho1[w]*G[w][w];
879 temp_PDF2 += t1D_2[w]*t1rho2[w]*G[w][w];
880 }
881 tmp1 = temp_PDF1*B1_rho1*B1_D_Kstr2_rho1;
882 tmp2 = temp_PDF2*B1_rho2*B1_D_Kstr1_rho2;
883 if(g1[i]==1) {
884 propagatorGS(mass2sq,mass2[i],width2[i],srho1,spion1,spi0,rRes2,pro_1);
885 propagatorGS(mass2sq,mass2[i],width2[i],srho2,spion2,spi0,rRes2,pro_2);
886 } else if(g1[i]==0) {
887 pro_1[0] = 1; pro_1[1] = 0; pro_2[0] = 1; pro_2[1] = 0;
888 }
889 temp_1[0] = tmp1*pro_1[0];
890 temp_1[1] = tmp1*pro_1[1];
891 temp_2[0] = tmp2*pro_2[0];
892 temp_2[1] = tmp2*pro_2[1];
893 KPiSLASS(skstr1,skaon,spion1,Amp_KPiS1);
894 KPiSLASS(skstr2,skaon,spion2,Amp_KPiS2);
895 Com_Multi(temp_1,Amp_KPiS2,amp_tmp1);
896 Com_Multi(temp_2,Amp_KPiS1,amp_tmp2);
897 }
898 else if (modetype[i]==9) {
899 if (B1_kstr1_Pi1<0.0) B1_kstr1_Pi1 = Barrier(1,skstr1,skaon,spion1,rRes2,mass1[i]);
900 if (B1_kstr2_Pi2<0.0) B1_kstr2_Pi2 = Barrier(1,skstr2,skaon,spion2,rRes2,mass1[i]);
901 if (B1_D_K1p_Pi0<0.0) B1_D_K1p_Pi0 = Barrier(1,sD,sk13,spi0,rD2,mD);
902 calt1(pK13,Pi0,t1D_1);
903 if(g2[i]==0) {
904 for(int w=0; w<4; w++) {
905 for(int j=0; j<4; j++) {
906 temp_PDF1 += t1D_1[w]*(G[w][j]-pK13[w]*pK13[j]/sk13)*t1kstr1[j]*G[w][w]*G[j][j];
907 temp_PDF2 += t1D_1[w]*(G[w][j]-pK13[w]*pK13[j]/sk13)*t1kstr2[j]*G[w][w]*G[j][j];
908 }
909 }
910 tmp1 = B1_kstr1_Pi1*B1_D_K1p_Pi0*temp_PDF1;
911 tmp2 = B1_kstr2_Pi2*B1_D_K1p_Pi0*temp_PDF2;
912 } else if(g2[i]==2) {
913 for(int w=0; w<4; w++) {
914 for(int j=0; j<4; j++) {
915 temp_PDF1 += t1D_1[w]*t2kk1[w][j]*t1kstr1[j]*G[w][w]*G[j][j];
916 temp_PDF2 += t1D_1[w]*t2kk2[w][j]*t1kstr2[j]*G[w][w]*G[j][j];
917 }
918 }
919 if (B2_K1p_kstr1<0.0) B2_K1p_kstr1 = Barrier(2,sk13,skstr1,spion2,rRes2,mass2[i]);
920 if (B2_K1p_kstr2<0.0) B2_K1p_kstr2 = Barrier(2,sk13,skstr2,spion1,rRes2,mass2[i]);
921 tmp1 = B1_kstr1_Pi1*B2_K1p_kstr1*B1_D_K1p_Pi0*temp_PDF1;
922 tmp2 = B1_kstr2_Pi2*B2_K1p_kstr2*B1_D_K1p_Pi0*temp_PDF2;
923 }
924 if(g0[i]==1) {
925 propagatorRBW(mass1sq,mass1[i],width1[i],skstr1,skaon,spion1,rRes2,1,pro0_1);
926 propagatorRBW(mass1sq,mass1[i],width1[i],skstr2,skaon,spion2,rRes2,1,pro0_2);
927 } else if(g0[i]==0) {
928 pro0_1[0] = 1; pro0_1[1] = 0; pro0_2[0] = 1; pro0_2[1] = 0;
929 }
930 if(g1[i]==1) {
931 propagatorRBW(mass2sq,mass2[i],width2[i],sk13,skstr1,spion2,rRes2,g2[i],pro1_1);
932 propagatorRBW(mass2sq,mass2[i],width2[i],sk13,skstr2,spion1,rRes2,g2[i],pro1_2);
933 } else if(g1[i]==0) {
934 pro1_1[0] = 1; pro1_1[1] = 0; pro1_2[0] = 1; pro1_2[1] = 0;
935 }
936 Com_Multi(pro0_1,pro1_1,pro_1);
937 Com_Multi(pro0_2,pro1_2,pro_2);
938 amp_tmp1[0] = tmp1*pro_1[0];
939 amp_tmp1[1] = tmp1*pro_1[1];
940 amp_tmp2[0] = tmp2*pro_2[0];
941 amp_tmp2[1] = tmp2*pro_2[1];
942 }
943
944 amp_tmp[0] = amp_tmp1[0]+amp_tmp2[0];
945 amp_tmp[1] = amp_tmp1[1]+amp_tmp2[1];
946
947 Com_Multi(amp_tmp,cof,amp_PDF);
948 PDF[0] += amp_PDF[0];
949 PDF[1] += amp_PDF[1];
950 double value = PDF[0]*PDF[0] + PDF[1]*PDF[1];
951 if(value <=0) {value = 1e-20;}
952 Result = value;
953
954}
955}
sin
double sin(const BesAngle a)
Definition: BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: BesAngle.h:213
mass
double mass
Definition: CosmicGenerator.cxx:138
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
g1
TF1 * g1
Definition: DataBase/tau_mode.c:66
w
double w
Definition: EvtPhokharaDef.hh:29
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
EvtDToKmPipPipPi0::decay
void decay(EvtParticle *p)
Definition: EvtDToKmPipPipPi0.cc:148
EvtDToKmPipPipPi0::getName
void getName(std::string &name)
Definition: EvtDToKmPipPipPi0.cc:39
EvtDToKmPipPipPi0::~EvtDToKmPipPipPi0
virtual ~EvtDToKmPipPipPi0()
Definition: EvtDToKmPipPipPi0.cc:37
EvtDToKmPipPipPi0::clone
EvtDecayBase * clone()
Definition: EvtDToKmPipPipPi0.cc:43
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
EvtPDL::getStdHep
static int getStdHep(EvtId id)
Definition: EvtPDL.hh:56
EvtParticle::getId
EvtId getId() const
Definition: EvtParticle.cc:112
EvtParticle::getP4
const EvtVector4R & getP4() const
Definition: EvtParticle.cc:120
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition: EvtParticle.cc:84
EvtParticle::initializePhaseSpace
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
Definition: EvtParticle.cc:1070
EvtVector4R::get
double get(int i) const
Definition: EvtVector4R.hh:179
b
const double b
Definition: slope.cxx:9