BOSS 7.1.2
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDsToKpPipPimPi0.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2// Environment:
3// This software is part of the EvtGen package developed jointly
4// for the BaBar and CLEO collaborations. If you use all or part
5// of it, please give an appropriate acknowledgement.
6//
7// Copyright Information: See EvtGen/COPYRIGHT
8// Copyright (C) 1998 Caltech, UCSB
9//
10// Module: EvtDsToKpPipPimPi0.cc
11//
12// Description: Routine to handle three-body decays of Ds+/Ds-
13//
14// Modification history:
15// Liaoyuan Dong Sun Jun 02 01:23:25 2024 Module created
16//------------------------------------------------------------------------
17#include "EvtGenBase/EvtPatches.hh"
18#include <fstream>
19#include <stdlib.h>
20#include "EvtGenBase/EvtParticle.hh"
21#include "EvtGenBase/EvtGenKine.hh"
22#include "EvtGenBase/EvtPDL.hh"
23#include "EvtGenBase/EvtReport.hh"
24#include "EvtGenBase/EvtResonance.hh"
25#include "EvtGenBase/EvtResonance2.hh"
26#include "EvtGenModels/EvtDsToKpPipPimPi0.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
37EvtDsToKpPipPimPi0::~EvtDsToKpPipPimPi0() {}
38
39void EvtDsToKpPipPimPi0::getName(std::string& model_name){
40 model_name="DsToKpPipPimPi0";
41}
42
46
47void EvtDsToKpPipPimPi0::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 int mode = 0;
59 phi[1] = -4.1936e+00;
60 phi[2] = 2.4056e+00;
61 phi[3] = 1.8102e+00;
62 phi[4] = -1.6106e+00;
63 phi[5] = -1.6106e+00;
64 phi[6] = 5.1214e+00;
65 phi[7] = 5.1214e+00;
66 phi[8] = 6.0675e-01;
67 phi[9] = -5.2450e+00;
68 phi[10]= -2.8735e+00;
69
70 rho[1] = 1.0981e+00;
71 rho[2] = 4.6322e-01;
72 rho[3] = 1.3138e+00;
73 rho[4] = 1.2709e+00;
74 rho[5] = 1.2709e+00;
75 rho[6] = 1.2953e+00;
76 rho[7] = 1.2953e+00;
77 rho[8] = 2.3583e+00;
78 rho[9] = 7.8907e+00;
79 rho[10]= 6.5982e-01;
80
81 phi[0] = 0;
82 rho[0] = 1;
83 mDsM = 1.9683;
84
85//------------------------new----------------------
86 mKst0 = 0.89555;
87 mKstp = 0.89166;
88 mrho = 0.7751;
89 mrho0 = 0.7753;
90 mK1270 = 1.289;
91 mK1400 = 1.403;
92 mK1650 = 1.672;
93 mOmega = 0.78266;
94 mA1 = 1.230;
95
96 GKst0 = 0.0473;
97 GKstp = 0.0508;
98 Grho = 0.1491;
99 Grho0 = 0.1478;
100 GK1270 = 0.116;
101 GK1400 = 0.174;
102 GK1650 = 0.158;
103 GOmega = 0.00868;
104 GA1 = 0.42;
105
106 //cout << "DsToKpPipPimPi0 :" << endl;
107 //for (int i=0; i<10; i++) {
108 // cout << i << " rho= " << rho[i] << " phi= " << phi[i] << endl;
109 //}
110
111 mass_Pion = 0.13957;
112 mass_Pion_N = 0.134977;
113 mass_Eta = 0.547862;
114 mass_Kaon = 0.493677;
115 math_pi = 3.1415926;
116
117 rD2 = 5.0; // 5*5
118 rRes1 = 3.0; // 3*3
119 rRes2 = 9.0; // 3*3
120
121 GS1 = 0.636619783;
122 GS2 = 0.01860182466;
123 GS3 = 0.1591549458; // 1/(2*math_2pi)
124 GS4 = 0.00620060822; // mass_Pion2/math_pi
125
126 int GG[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
127 int EE[4][4][4][4] =
128 { { {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
129 {{0,0,0,0}, {0,0,0,0}, {0,0,0,1}, {0,0,-1,0}},
130 {{0,0,0,0}, {0,0,0,-1}, {0,0,0,0}, {0,1,0,0} },
131 {{0,0,0,0}, {0,0,1,0}, {0,-1,0,0}, {0,0,0,0} } },
132 { {{0,0,0,0}, {0,0,0,0}, {0,0,0,-1}, {0,0,1,0} },
133 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
134 {{0,0,0,1}, {0,0,0,0}, {0,0,0,0}, {-1,0,0,0}},
135 {{0,0,-1,0}, {0,0,0,0}, {1,0,0,0}, {0,0,0,0} } },
136 { {{0,0,0,0}, {0,0,0,1}, {0,0,0,0}, {0,-1,0,0}},
137 {{0,0,0,-1}, {0,0,0,0}, {0,0,0,0}, {1,0,0,0} },
138 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
139 {{0,1,0,0}, {-1,0,0,0}, {0,0,0,0}, {0,0,0,0} } },
140 { {{0,0,0,0}, {0,0,-1,0}, {0,1,0,0}, {0,0,0,0} },
141 {{0,0,1,0}, {0,0,0,0}, {-1,0,0,0}, {0,0,0,0} },
142 {{0,-1,0,0}, {1,0,0,0}, {0,0,0,0}, {0,0,0,0} },
143 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} } } };
144 for (int i=0; i<4; i++) {
145 for (int j=0; j<4; j++) {
146 G[i][j] = GG[i][j];
147 for (int k=0; k<4; k++) {
148 for (int l=0; l<4; l++) {
149 E[i][j][k][l] = EE[i][j][k][l];
150 }
151 }
152 }
153 }
154}
155
156
160
161void EvtDsToKpPipPimPi0::decay( EvtParticle *p){
162//--------------max
163
164/* double maxprob=0,maxq1270;
165 for(int ir=0;ir<=60000000;ir++){
166 p->initializePhaseSpace(getNDaug(),getDaugs());
167 EvtVector4R _kp = p->getDaug(0)->getP4();
168 EvtVector4R _pip = p->getDaug(1)->getP4();
169 EvtVector4R _pim = p->getDaug(2)->getP4();
170 EvtVector4R _pi0 = p->getDaug(3)->getP4();
171
172 double _Pip[4],_Pim[4],_Pi0[4],_Kp[4];
173 _Kp[0] = _kp.get(0); _Pip[0] = _pip.get(0); _Pim[0] = _pim.get(0); _Pi0[0] = _pi0.get(0);
174 _Kp[1] = _kp.get(1); _Pip[1] = _pip.get(1); _Pim[1] = _pim.get(1); _Pi0[1] = _pi0.get(1);
175 _Kp[2] = _kp.get(2); _Pip[2] = _pip.get(2); _Pim[2] = _pim.get(2); _Pi0[2] = _pi0.get(2);
176 _Kp[3] = _kp.get(3); _Pip[3] = _pip.get(3); _Pim[3] = _pim.get(3); _Pi0[3] = _pi0.get(3);
177
178 double _prob;
179 double value,q1270;
180 int nstates=11;
181 int modetype[11]= {1,1,2,10,12,11,16,17,14,2,15};
182 int g0[11] = {1,1,1,1,1,1,1,1,1,0,0};
183 int g1[11] = {1,1,1,1,1,1,1,1,1,1,0};
184 int g2[11] = {0,1,1,0,0,0,0,0,1,0,0};
185 double mass1[11] = {mKst0,mKst0,mKstp,mK1270,mKstp,mKst0,mA1,mA1,mOmega,mKst0,mKst0};
186 double mass2[11] = {mrho,mrho,mrho0,mrho,mK1400,mK1400,mrho,mrho,mrho,mrho,mrho};
187 double width1[11] = {GKst0,GKst0,GKstp,GK1270,GKstp,GKst0,GA1,GA1,GOmega,GKst0,GKst0};
188 double width2[11] = {Grho,Grho,Grho0,Grho,GK1400,GK1400,Grho,Grho,Grho,mrho,mrho};
189
190 // calEvaMy(_Kp,_Pip,_Pim,_Pi0,mass1,mass2,width1,width2,rho,phi,g0,g1,g2,modetype,nstates,_prob);
191 calEvaMy(_Kp,_Pip,_Pim,_Pi0,mass1,mass2,width1,width2,rho,phi,g0,g1,g2,modetype,nstates,_prob,q1270);
192
193 if(_prob>maxprob) {
194 maxprob=_prob;
195 maxq1270=q1270;
196 printf("ir = %d,maxprob = %.10f,prob = %.10f,q1270 = %.10f\n\n",ir,maxprob,_prob,maxq1270);
197// cout<<"mm: "<<_Kp[0]<<" "<<_Pip[0]<<" "<<_Pim[0]<<" "<<_Pi0[0]<<endl;
198// cout<<"mm: "<<_Kp[1]<<" "<<_Pip[1]<<" "<<_Pim[1]<<" "<<_Pi0[1]<<endl;
199// cout<<"mm: "<<_Kp[2]<<" "<<_Pip[2]<<" "<<_Pim[2]<<" "<<_Pi0[2]<<endl;
200// cout<<"mm: "<<_Kp[3]<<" "<<_Pip[3]<<" "<<_Pim[3]<<" "<<_Pi0[3]<<endl;
201 }
202 }
203 printf("maxprob = %.10f\n", maxprob);
204*/
205
206 p->initializePhaseSpace(getNDaug(),getDaugs());
207 EvtVector4R kp = p->getDaug(0)->getP4();
208 EvtVector4R pip = p->getDaug(1)->getP4();
209 EvtVector4R pim = p->getDaug(2)->getP4();
210 EvtVector4R pi0 = p->getDaug(3)->getP4();
211
212 double Pip[4],Pim[4],Kp[4],Pi0[4];
213
214 //----------------------------------------------------------------------------------------
215 Kp[0] = kp.get(0); Pip[0] = pip.get(0); Pim[0] = pim.get(0);Pi0[0] = pi0.get(0);
216 Kp[1] = kp.get(1); Pip[1] = pip.get(1); Pim[1] = pim.get(1);Pi0[1] = pi0.get(1);
217 Kp[2] = kp.get(2); Pip[2] = pip.get(2); Pim[2] = pim.get(2);Pi0[2] = pi0.get(2);
218 Kp[3] = kp.get(3); Pip[3] = pip.get(3); Pim[3] = pim.get(3);Pi0[3] = pi0.get(3);
219 //---------------------------------------------------------------------------------------
220 double value,q1270;
221 int nstates=11;
222 int modetype[11]= {1,1,2,10,12,11,16,17,14,2,15};
223 int g0[11] = {1,1,1,1,1,1,1,1,1,0,0};
224 int g1[11] = {1,1,1,1,1,1,1,1,1,1,0};
225 int g2[11] = {0,1,1,0,0,0,0,0,1,0,0};
226 double mass1[11] = {mKst0,mKst0,mKstp,mK1270,mKstp,mKst0,mA1,mA1,mOmega,mKst0,mKst0};
227 double mass2[11] = {mrho,mrho,mrho0,mrho,mK1400,mK1400,mrho,mrho,mrho,mrho,mrho};
228 double width1[11] = {GKst0,GKst0,GKstp,GK1270,GKstp,GKst0,GA1,GA1,GOmega,GKst0,GKst0};
229 double width2[11] = {Grho,Grho,Grho0,Grho,GK1400,GK1400,Grho,Grho,Grho,mrho,mrho};
230
231 calEvaMy(Kp,Pip,Pim,Pi0,mass1,mass2,width1,width2,rho,phi,g0,g1,g2,modetype,nstates,value,q1270);
232 setProb(value);
233 return;
234}
235
236double EvtDsToKpPipPimPi0::CalRho4pi(double s)
237 {
238 if(s>=1.){
239 return sqrt((s-16.*mass_Pion*mass_Pion)/s);
240 }else{
241 double gam = 0.0005;
242 double s2 = s*s;
243 double s3 = s2*s;
244 double s4 = s2*s2;
245 double s5 = s4*s;
246 double s6 = s5*s;
247
248 gam -= 0.0193*s;
249 gam += 0.1385*s2;
250 gam -= 0.2084*s3;
251 gam -= 0.2974*s4;
252 gam += 0.1366*s5;
253 gam += 1.0789*s6;
254
255 return gam;
256 }
257 }
258TComplex EvtDsToKpPipPimPi0::ResonanceSkm(double & m2)
259{
260 double g11 = 0.22889, g12 = -0.55377, g13 = 0.00, g14 = -0.39899, g15 = -0.34639;
261 double g21 = 0.94128, g22 = 0.55095, g23 = 0.00, g24 = 0.39065, g25 = 0.31503;
262 double g31 = 0.36856, g32 = 0.23888, g33 = 0.55639, g34 = 0.18340, g35 = 0.18681;
263 double g41 = 0.33650, g42 = 0.40907, g43 = 0.85679, g44 = 0.19906, g45 = -0.00984;
264 double g51 = 0.18171, g52 = -0.17558, g53 = -0.79658, g54 = -0.00355, g55 = 0.22358;
265 double fs11 = 0.23399, fs12 = 0.15044, fs13 =-0.20545, fs14 = 0.32825, fs15 = 0.35412;
266 double m_1 = 0.65100, m_2 = 1.20360, m_3 = 1.55817, m_4 = 1.21000, m_5 = 1.82206;
267 double ss0 = -3.92637;
268 double sp0 = -3.0;//v0
269 double sA0 = -0.15;
270 double sA = 1.0;
271 const double mass_Eta = 0.547862;
272 const double mass_Etap = 0.95778;
273
274 double km11 = (g11*g11/(m_1*m_1-m2)+g21*g21/(m_2*m_2-m2)+g31*g31/(m_3*m_3-m2)+g41*g41/(m_4*m_4-m2)+g51*g51/(m_5*m_5-m2)+fs11*(1-ss0)/(m2-ss0))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
275 double km12 = (g11*g12/(m_1*m_1-m2)+g21*g22/(m_2*m_2-m2)+g31*g32/(m_3*m_3-m2)+g41*g42/(m_4*m_4-m2)+g51*g52/(m_5*m_5-m2)+fs12*(1-ss0)/(m2-ss0))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
276 double km13 = (g11*g13/(m_1*m_1-m2)+g21*g23/(m_2*m_2-m2)+g31*g33/(m_3*m_3-m2)+g41*g43/(m_4*m_4-m2)+g51*g53/(m_5*m_5-m2)+fs13*(1-ss0)/(m2-ss0))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
277 double km14 = (g11*g14/(m_1*m_1-m2)+g21*g24/(m_2*m_2-m2)+g31*g34/(m_3*m_3-m2)+g41*g44/(m_4*m_4-m2)+g51*g54/(m_5*m_5-m2)+fs14*(1-ss0)/(m2-ss0))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
278 double km15 = (g11*g15/(m_1*m_1-m2)+g21*g25/(m_2*m_2-m2)+g31*g35/(m_3*m_3-m2)+g41*g45/(m_4*m_4-m2)+g51*g55/(m_5*m_5-m2)+fs15*(1-ss0)/(m2-ss0))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
279 double km21 = km12, km31 = km13, km41 = km14, km51 = km15;
280 double km23 = (g12*g13/(m_1*m_1-m2)+g22*g23/(m_2*m_2-m2)+g32*g33/(m_3*m_3-m2)+g42*g43/(m_4*m_4-m2)+g52*g53/(m_5*m_5-m2))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
281 double km24 = (g12*g14/(m_1*m_1-m2)+g22*g24/(m_2*m_2-m2)+g32*g34/(m_3*m_3-m2)+g42*g44/(m_4*m_4-m2)+g52*g54/(m_5*m_5-m2))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
282 double km25 = (g12*g15/(m_1*m_1-m2)+g22*g25/(m_2*m_2-m2)+g32*g35/(m_3*m_3-m2)+g42*g45/(m_4*m_4-m2)+g52*g55/(m_5*m_5-m2))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
283 double km32 = km23, km42 = km24, km52 = km25;
284 double km34 = (g13*g14/(m_1*m_1-m2)+g23*g24/(m_2*m_2-m2)+g33*g34/(m_3*m_3-m2)+g43*g44/(m_4*m_4-m2)+g53*g54/(m_5*m_5-m2))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
285 double km35 = (g13*g15/(m_1*m_1-m2)+g23*g25/(m_2*m_2-m2)+g33*g35/(m_3*m_3-m2)+g43*g45/(m_4*m_4-m2)+g53*g55/(m_5*m_5-m2))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
286 double km43 = km34, km53 = km35;
287 double km45 = (g14*g15/(m_1*m_1-m2)+g24*g25/(m_2*m_2-m2)+g34*g35/(m_3*m_3-m2)+g44*g45/(m_4*m_4-m2)+g54*g55/(m_5*m_5-m2))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
288 double km54 = km45;
289 double km22 = (g12*g12/(m_1*m_1-m2)+g22*g22/(m_2*m_2-m2)+g32*g32/(m_3*m_3-m2)+g42*g42/(m_4*m_4-m2)+g52*g52/(m_5*m_5-m2))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
290 double km33 = (g13*g13/(m_1*m_1-m2)+g23*g23/(m_2*m_2-m2)+g33*g33/(m_3*m_3-m2)+g43*g43/(m_4*m_4-m2)+g53*g53/(m_5*m_5-m2))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
291 double km44 = (g14*g14/(m_1*m_1-m2)+g24*g24/(m_2*m_2-m2)+g34*g34/(m_3*m_3-m2)+g44*g44/(m_4*m_4-m2)+g54*g54/(m_5*m_5-m2 ))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
292 double km55 = (g15*g15/(m_1*m_1-m2)+g25*g25/(m_2*m_2-m2)+g35*g35/(m_3*m_3-m2)+g45*g45/(m_4*m_4-m2)+g55*g55/(m_5*m_5-m2 ))*(1-sA0)/(m2-sA0)*(m2-sA*mass_Pion*mass_Pion*0.5);
293 TComplex fp11(16.55, -3.87); TComplex beta1(8.89, 3.36);
294 TComplex fp12(-13.99, 5.12); TComplex beta2(16.20, 5.77);
295 TComplex fp13(-77.34, -80.04); TComplex beta3(-21.60, 27.41);
296 TComplex fp14(-28.52, -3.20); TComplex beta4(-40.15, 35.36);
297 TComplex fp15(30.47, 5.70); TComplex beta5(30.0, -43.09);
298 TComplex P1 = fp11*(1-sp0)/(m2-sp0)+beta1*g11/(m_1*m_1-m2)+beta2*g21/(m_2*m_2-m2)+beta3*g31/(m_3*m_3-m2)+beta4*g41/(m_4*m_4-m2)+beta5*g51/(m_5*m_5-m2);
299 TComplex P2 = fp12*(1-sp0)/(m2-sp0)+beta1*g12/(m_1*m_1-m2)+beta2*g22/(m_2*m_2-m2)+beta3*g32/(m_3*m_3-m2)+beta4*g42/(m_4*m_4-m2)+beta5*g52/(m_5*m_5-m2);
300 TComplex P3 = fp13*(1-sp0)/(m2-sp0)+beta1*g13/(m_1*m_1-m2)+beta2*g23/(m_2*m_2-m2)+beta3*g33/(m_3*m_3-m2)+beta4*g43/(m_4*m_4-m2)+beta5*g53/(m_5*m_5-m2);
301 TComplex P4 = fp14*(1-sp0)/(m2-sp0)+beta1*g14/(m_1*m_1-m2)+beta2*g24/(m_2*m_2-m2)+beta3*g34/(m_3*m_3-m2)+beta4*g44/(m_4*m_4-m2)+beta5*g54/(m_5*m_5-m2);
302 TComplex P5 = fp15*(1-sp0)/(m2-sp0)+beta1*g15/(m_1*m_1-m2)+beta2*g25/(m_2*m_2-m2)+beta3*g35/(m_3*m_3-m2)+beta4*g45/(m_4*m_4-m2)+beta5*g55/(m_5*m_5-m2);
303
304 TMatrixD MI(5,5), MA(5,5), MA_invt(5,5), MB(5,5), KM(5,5), RhoA(5,5), RhoB(5,5), MRe(5,5), MIm(5,5);
305 MI.UnitMatrix();
306
307 RhoA.UnitMatrix();
308 RhoB.UnitMatrix();
309 TMatrixDRow(KM,0)(0) = km11;TMatrixDRow(KM,1)(0) = km21;TMatrixDRow(KM,2)(0) = km31;TMatrixDRow(KM,3)(0) = km41;TMatrixDRow(KM,4)(0)= km51;
310 TMatrixDRow(KM,0)(1) = km12;TMatrixDRow(KM,1)(1) = km22;TMatrixDRow(KM,2)(1) = km32;TMatrixDRow(KM,3)(1) = km42;TMatrixDRow(KM,4)(1)= km52;
311 TMatrixDRow(KM,0)(2) = km13;TMatrixDRow(KM,1)(2) = km23;TMatrixDRow(KM,2)(2) = km33;TMatrixDRow(KM,3)(2) = km43;TMatrixDRow(KM,4)(2)= km53;
312 TMatrixDRow(KM,0)(3) = km14;TMatrixDRow(KM,1)(3) = km24;TMatrixDRow(KM,2)(3) = km34;TMatrixDRow(KM,3)(3) = km44;TMatrixDRow(KM,4)(3)= km54;
313 TMatrixDRow(KM,0)(4) = km15;TMatrixDRow(KM,1)(4) = km25;TMatrixDRow(KM,2)(4) = km35;TMatrixDRow(KM,3)(4) = km45;TMatrixDRow(KM,4)(4)= km55;
314
315 TMatrixDRow(RhoA,0)(0) = sqrt(1.0-4.0*mass_Pion*mass_Pion/m2);
316 TMatrixDRow(RhoB,0)(0) = 0.0;
317 if ( (1.0-4.0*mass_Kaon*mass_Kaon/m2) > 0) {
318 TMatrixDRow(RhoA,1)(1) = sqrt(1.0-4.0*mass_Kaon*mass_Kaon/m2);
319 TMatrixDRow(RhoB,1)(1) = 0.0;
320 }
321 else {
322 TMatrixDRow(RhoA,1)(1) = 0.0;
323 TMatrixDRow(RhoB,1)(1) = sqrt(4.0*mass_Kaon*mass_Kaon/m2-1.0);
324 }
325 TMatrixDRow(RhoA,2)(2) = CalRho4pi(m2);
326 TMatrixDRow(RhoB,2)(2) = 0.0;
327 if ( (1.0-4.0*mass_Eta*mass_Eta/m2) > 0) {
328 TMatrixDRow(RhoA,3)(3) = sqrt(1.0-4.0*mass_Eta*mass_Eta/m2);
329 TMatrixDRow(RhoB,3)(3) = 0.0;
330 }else {
331 TMatrixDRow(RhoA,3)(3) = 0.0;
332 TMatrixDRow(RhoB,3)(3) = sqrt(4.0*mass_Eta*mass_Eta/m2-1.0);
333 }
334 TMatrixDRow(RhoA,4)(4) = 0.0;
335 TMatrixDRow(RhoB,4)(4) = sqrt((mass_Eta + mass_Etap)*(mass_Eta + mass_Etap)/m2-1.0);
336
337 MA = MI + KM*RhoB;
338 MB = -1.0*KM*RhoA;
339 MA_invt = MA;
340 MA_invt.Invert();
341 MRe = MA+MB*MA_invt*MB;
342 MRe.Invert();
343 MIm = MA_invt*MB*MRe;
344 TComplex ciR(1.0, 0.0);
345 TComplex ciM(0.0, 1.0);
346 TComplex amp;
347 amp = (ciR*TMatrixDRow(MRe,0)(0)-ciM*TMatrixDRow(MIm,0)(0))*P1+
348 (ciR*TMatrixDRow(MRe,0)(1)-ciM*TMatrixDRow(MIm,0)(1))*P2+
349 (ciR*TMatrixDRow(MRe,0)(2)-ciM*TMatrixDRow(MIm,0)(2))*P3+
350 (ciR*TMatrixDRow(MRe,0)(3)-ciM*TMatrixDRow(MIm,0)(3))*P4+
351 (ciR*TMatrixDRow(MRe,0)(4)-ciM*TMatrixDRow(MIm,0)(4))*P5;
352 return amp;
353 }
354void EvtDsToKpPipPimPi0::KPiSLASS(double sa, double sb, double sc, double prop[2]) {
355 const double m1430 = 1.463;
356 const double sa0 = 2.140369; // m1430*m1430;
357 const double w1430 = 0.233;
358 const double Lass1 = 0.25/sa0;
359 double tmp = sb-sc;
360 double tmp1 = sa0+tmp;
361 double q0 = Lass1*tmp1*tmp1-sb;
362 if(q0<0) q0 = 1e-16;
363 double tmp2 = sa+tmp;
364 double qs = 0.25*tmp2*tmp2/sa-sb;
365 double q = sqrt(qs);
366 double width = w1430*q*m1430/sqrt(sa*q0);
367 double temp_R = atan(m1430*width/(sa0-sa));
368 if(temp_R<0) temp_R += math_pi;
369 double deltaR = -5.31 + temp_R;
370 double temp_F = atan(2.14*q/(2.0-1.926*qs)); // 2.0*1.07 = 2.14; 1.8*1.07 = 1.926
371 if(temp_F<0) temp_F += math_pi;
372 double deltaF = 2.33 + temp_F;
373 double deltaS = deltaR + 2.0*deltaF;
374 double t1 = 0.8*sin(deltaF);
375 double t2 = sin(deltaR);
376 double CF[2], CS[2];
377 CF[0] = cos(deltaF);
378 CF[1] = sin(deltaF);
379 CS[0] = cos(deltaS);
380 CS[1] = sin(deltaS);
381 prop[0] = t1*CF[0] + t2*CS[0];
382 prop[1] = t1*CF[1] + t2*CS[1];
383}
384
385
386void EvtDsToKpPipPimPi0::Com_Multi(double a1[2], double a2[2], double res[2]){
387 res[0] = a1[0]*a2[0]-a1[1]*a2[1];
388 res[1] = a1[1]*a2[0]+a1[0]*a2[1];
389}
390void EvtDsToKpPipPimPi0::Com_Divide(double a1[2], double a2[2], double res[2]){
391 res[0] = (a1[0]*a2[0]+a1[1]*a2[1])/(a2[0]*a2[0]+a2[1]*a2[1]);
392 res[1] = (a1[1]*a2[0]-a1[0]*a2[1])/(a2[0]*a2[0]+a2[1]*a2[1]);
393}
394double EvtDsToKpPipPimPi0::SCADot(double a1[4], double a2[4])
395{
396 double _cal = 0.;
397 _cal = a1[0]*a2[0]-a1[1]*a2[1]-a1[2]*a2[2]-a1[3]*a2[3];
398 return _cal;
399}
400double EvtDsToKpPipPimPi0::barrier(int l, double sa, double sb, double sc, double r,double mass)
401{
402 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
403// if(q < 0) q = 1e-16;
404 if(q < 0) q = -q;
405
406 double z;
407 z = q*r*r;
408 double sa0;
409 sa0 = mass*mass;
410 double q0 = (sa0+sb-sc)*(sa0+sb-sc)/(4*sa0)-sb;
411// if(q0 < 0) q0 = 1e-16;
412 if(q0 < 0) q0 = -q0;
413 double z0 = q0*r*r;
414 double F = 0.0;
415 if(l == 0) F = 1;
416 if(l == 1) F = sqrt((1+z0)/(1+z));
417 if(l == 2) F = sqrt((9+3*z0+z0*z0)/(9+3*z+z*z));
418 // if(l == 1) F = sqrt((2*r)/(1+z));
419 // if(l == 2) F = sqrt((13*r*r)/(9+3*z+z*z));
420 return F;
421}
422void EvtDsToKpPipPimPi0::calt1(double daug1[4], double daug0[4], double t1[4]){
423 double p, pq;
424 double pa[4], qa[4];
425 for(int i=0; i<4; i++){
426 pa[i] = daug1[i] + daug0[i];
427 qa[i] = daug1[i] - daug0[i];
428 }
429 p = SCADot(pa,pa);
430 pq = SCADot(pa,qa);
431 for(int i=0; i<4; i++){
432 t1[i] = qa[i] - pq/p*pa[i];
433 }
434}
435void EvtDsToKpPipPimPi0::calt2(double daug1[4], double daug0[4], double t2[4][4]){
436 double p, r;
437 double pa[4], t1[4];
438 calt1(daug1,daug0,t1);
439 r = SCADot(t1,t1)/3.0;
440 for(int i=0; i<4; i++) {
441 pa[i] = daug1[i] + daug0[i];
442 }
443 p = SCADot(pa,pa);
444 for(int i=0; i<4; i++) {
445 for(int j=0; j<4; j++) {
446 t2[i][j] = t1[i]*t1[j] - r*(G[i][j]-pa[i]*pa[j]/p);
447 }
448 }
449}
450
451double EvtDsToKpPipPimPi0::wid(double mass2, double mass, double sa, double sb, double sc, double r, int l){
452 double widm = 0.;
453 double m = sqrt(sa);
454 double tmp = sb-sc;
455 double tmp1 = sa+tmp;
456 double q = 0.25*tmp1*tmp1/sa-sb;
457// if(q<0) q = 1e-16;
458 if(q<0) q = -q;
459 double tmp2 = mass2+tmp;
460 double q0 = 0.25*tmp2*tmp2/mass2-sb;
461// if(q0<0) q0 = 1e-16;
462 if(q0<0) q0 = -q0;
463 double z = q*r*r;
464 double z0 = q0*r*r;
465 double t = q/q0;
466 if(l == 0) {widm = sqrt(t)*mass/m;}
467 else if(l == 1) {widm = t*sqrt(t)*mass/m*(1+z0)/(1+z);}
468 else if(l == 2) {widm = t*t*sqrt(t)*mass/m*(9+3*z0+z0*z0)/(9+3*z+z*z);}
469 return widm;
470
471}
472double EvtDsToKpPipPimPi0::widl1(double mass2, double mass, double sa, double sb, double sc, double r)
473{
474 double widm = 0.;
475 double m = sqrt(sa);
476 double tmp = sb-sc;
477 double tmp1 = sa+tmp;
478 double q = 0.25*tmp1*tmp1/sa-sb;
479// if(q<0) q = 1e-16;
480 if(q<0) q = -q;
481 double tmp2 = mass2+tmp;
482 double q0 = 0.25*tmp2*tmp2/mass2-sb;
483// if(q0<0) q0 = 1e-16;
484 if(q0<0) q0 = -q0;
485 double z = q*r*r;
486 double z0 = q0*r*r;
487 double F = (1+z0)/(1+z);
488 double t = q/q0;
489 widm = t*sqrt(t)*mass/m*F;
490 return widm;
491
492}
493void EvtDsToKpPipPimPi0::propagatorNBW(double mass2, double mass, double width, double sa, double sb, double sc, double r, int l, double prop[2]){
494 double a[2], b[2];
495 a[0] = 1;
496 a[1] = 0;
497 b[0] = mass2-sa;
498 b[1] = -mass*width;
499 Com_Divide(a,b,prop);
500}
501
502void EvtDsToKpPipPimPi0::propagatorRBW(double mass2, double mass, double width, double sa, double sb, double sc, double r, int l, double prop[2]){
503 double a[2], b[2];
504 a[0] = 1;
505 a[1] = 0;
506 b[0] = mass2-sa;
507 b[1] = -mass*width*wid(mass2,mass,sa,sb,sc,r,l);
508 Com_Divide(a,b,prop);
509}
510void EvtDsToKpPipPimPi0::propagatorRBWl1(double mass2, double mass, double width, double sa, double sb, double sc, double r, double prop[2])
511 //firstly code for propagator_omg, now could use for propagatorRBW that l=1, and will be more quick.
512{
513 double a[2], b[2];
514 a[0] = 1;
515 a[1] = 0;
516 b[0] = mass2-sa;
517 b[1] = -mass*width*widl1(mass2,mass,sa,sb,sc,r);
518 Com_Divide(a,b,prop);
519}
520//------------GS---used by rho----------------------------
521void EvtDsToKpPipPimPi0::propagatorGS(double mass2, double mass, double width, double sa, double sb, double sc, double r, double prop[2])
522{
523 double a[2], b[2];
524 double tmp = sb-sc;
525 double tmp1 = sa+tmp;
526 double q2 = 0.25*tmp1*tmp1/sa-sb;
527// if(q2<0) q2 = 1e-16;
528 if(q2<0) q2 = -q2;
529
530 double tmp2 = mass2+tmp;
531 double q02 = 0.25*tmp2*tmp2/mass2-sb;
532// if(q02<0) q02 = 1e-16;
533 if(q02<0) q02 = -q02;
534
535 double q = sqrt(q2);
536 double q0 = sqrt(q02);
537 double m = sqrt(sa);
538 double q03 = q0*q02;
539 double tmp3 = log(mass+2*q0)+1.2760418309; // log(mass_2Pion) = 1.2760418309;
540
541 double h = GS1*q/m*(log(m+2*q)+1.2760418309);
542 double h0 = GS1*q0/mass*tmp3;
543 double dh = h0*(0.125/q02-0.5/mass2)+GS3/mass2;
544 double d = GS2/q02*tmp3+GS3*mass/q0-GS4*mass/q03;
545 double f = mass2/q03*(q2*(h-h0)+(mass2-sa)*q02*dh);
546
547 a[0] = 1.0+d*width/mass;
548 a[1] = 0.0;
549 b[0] = mass2-sa+width*f;
550 b[1] = -mass*width*widl1(mass2,mass,sa,sb,sc,r);
551 Com_Divide(a,b,prop);
552}
553void EvtDsToKpPipPimPi0::calEvaMy(
554 double* Kp, double* Pip, double* Pim, double* Pi0,
555 double* mass1, double* mass2,
556 double* width1, double* width2,
557 double* amp, double* phase,
558 int* g0, int* g1, int* g2,
559 int* modetype, int nstates, double & Result, double & q1270)
560
561{
562
563 TComplex tmpswave;
564 double spipi(0.);
565
566 spipi = (Pip[0]+Pim[0])*(Pip[0]+Pim[0])-
567 (Pip[1]+Pim[1])*(Pip[1]+Pim[1])-
568 (Pip[2]+Pim[2])*(Pip[2]+Pim[2])-
569 (Pip[3]+Pim[3])*(Pip[3]+Pim[3]);
570 tmpswave = ResonanceSkm(spipi);
571 double realpipis = tmpswave.Re();
572 double imagpipis = tmpswave.Im();
573
574 double pKstr0[4],pKstrp[4],prhop[4],prhom[4],prho0[4],pK11[4],pK12[4],pK13[4],pA1[4],pD[4],pomega[4];
575 for(int i=0;i!=4;i++){
576 prhop[i] =Pip[i]+Pi0[i];
577 prhom[i] =Pim[i]+Pi0[i];
578 prho0[i] =Pim[i]+Pip[i];
579 pKstr0[i] =Kp[i]+Pim[i];
580 pKstrp[i] =Kp[i]+Pi0[i];
581 pK11[i] =pKstr0[i]+Pi0[i];
582 pK12[i] =Kp[i]+Pim[i];
583 pK13[i] =pKstr0[i]+Pip[i];
584 pA1[i] =prhop[i]+Pim[i];
585 pomega[i] =Pim[i]+Pip[i]+Pi0[i];
586 pD[i] =pKstr0[i]+prhop[i];
587 }
588 double spi0,spionp,spionm,skaon,srhop,srhom,srho0,sKst0,sKstp,sk11,sk12,sk13,sa1,sD,somega;
589
590 skaon = SCADot(Kp,Kp);
591 spionp = SCADot(Pip,Pip);
592 spionm = SCADot(Pim,Pim);
593 spi0 = SCADot(Pi0,Pi0);
594
595 srhop = SCADot(prhop,prhop);
596 srhom = SCADot(prhom,prhom);
597 srho0 = SCADot(prho0,prho0);
598 somega = SCADot(pomega,pomega);
599
600 sKst0 = SCADot(pKstr0,pKstr0);
601 sKstp = SCADot(pKstrp,pKstrp);
602 sk11 = SCADot(pK11,pK11);
603 sk12 = SCADot(pK12,pK12);
604 sk13 = SCADot(pK13,pK13);
605 sa1 = SCADot(pA1,pA1);
606 sD = SCADot(pD,pD);
607 //-------------------------------------------------------------------------------------------------
608 double t2A1[4][4],t2A2[4][4],t2A3[4][4],t1rhom[4],t1rhop[4],t1rho0[4],t1Kst0[4],t1Kstp[4],t2k11[4][4],t2k12[4][4];
609 double t1K1_KPi[4],t2k21[4][4],t2k22[4][4],t1A3[4];
610 calt1(Pip,Pi0,t1rhop);
611 calt1(Pim,Pi0,t1rhom);
612 calt1(Pim,Pip,t1rho0);
613
614 calt1(Kp,Pim,t1Kst0);
615 calt1(Kp,Pi0,t1Kstp);
616 calt1(Kp,Pim,t1K1_KPi);
617
618 calt1(prho0,Pi0,t1A3);
619
620 calt2(prhop,Pim,t2A1);
621 calt2(prhom,Pip,t2A2);
622 calt2(prho0,Pi0,t2A3);
623
624 calt2(pKstr0,Pi0,t2k11);
625 calt2(pKstrp,Pim,t2k12);
626
627 calt2(prhom,Kp,t2k21);
628 calt2(prho0,Kp,t2k22);
629
630 double B_rhop =-1.0, B_rhom =-1.0, B_Kst0rhop1 =-1.0, B_Kst0rhop2=-1.0;
631 double B_Kst0 =-1.0, B_K1_KPi =-1.0, B_K1rhop1 =-1.0, B_K1rhop2 =-1.0;
632 double B_Kstp =-1.0, B_K1rho =-1.0, B_D_K1rho =-1.0, B_RhopKs1 =-1.0;
633 double B_K1_Kst0pi0=-1.0, B_D_K1Pi =-1.0, B_K1_Kstppim=-1.0,B_K1Kstp2 =-1.0 ;
634 double B_D_A1 = -1.0, B_D_A2 =-1.0, B_rho0 =-1.0, B_Rho0Ks1 =-1.0;
635 double B_K1PiPi = -1.0, B_K14_1 =-1.0, B_K14_2 =-1.0;
636 double B_Kstprho01 =-1.0, B_Kstprho02=-1.0;
637
638 double cof[2],amp_tmp[2],amp_PDF[2], PDF[2];
639 PDF[0] = 0;
640 PDF[1] = 0;
641 double mass1sq,mass2sq,tt1,tt2,tt21,tt22,tt23,tt3,qqq=10,qqq0=10,nq=0;
642 double temp_PDF, tmp1,tmp3,temp_PDF1;
643 double pro[2], pro0[2], pro1[2],pro2[2],pro3[2],pro4[2];
644 double t1D[4],t2D[4][4],t1Ds_omega[4],B1_Ds_omega,B1_1V23;
645 int isKstp=0.0, isKst0=0.0, isRhop=0.0,isRhom=0.0,isRho0=0.0,isKpPim_S=0.0,isPipPi0_S=0.0,isf0=0.0,isPipPim_S=0.0;
646 double proRhop[2], proRhom[2],proRho0[2],proKstp[2], proKst0[2],proKPi_S[2],proPiPi_S[2],prof0[2];
647 double temp_PDF11, temp_PDF12, temp_PDF13;
648 double pro0_11[2],pro0_12[2],pro0_13[2];
649 double pro_11[2],pro_12[2],pro_13[2];
650
651 for(int i=0; i<nstates; i++){
652 temp_PDF = 0;
653 temp_PDF1 = 0;
654 temp_PDF11 = 0;
655 temp_PDF12 = 0;
656 temp_PDF13 = 0;
657 cof[0] = amp[i]*cos(phase[i]);
658 cof[1] = amp[i]*sin(phase[i]);
659 mass1sq = mass1[i]*mass1[i];//avoid parameters input reversed
660 mass2sq = mass2[i]*mass2[i];
661
662 if(modetype[i]==1){
663 if (B_Kst0<0.0) B_Kst0 = barrier(1,sKst0,skaon,spionm,rRes1,mass1[i]);
664 if (B_rhop<0.0) B_rhop = barrier(1,srhop,spionp,spi0,rRes1,mass2[i]);
665 if(g2[i]==0){
666 for(int w=0; w<4; w++){
667 temp_PDF += G[w][w]*t1Kst0[w]*t1rhop[w];
668 }
669 tmp1 = B_Kst0*B_rhop*temp_PDF;
670 }
671 else if(g2[i]==1){
672 calt1(pKstr0,prhop,t1D);
673 for(int w=0; w<4; w++){
674 tt1=pD[w]*G[w][w];
675 for(int j=0; j<4; j++){
676 tt2=t1D[j]*G[j][j];
677 for(int k=0; k<4; k++){
678 tt3=t1Kst0[k]*G[k][k];
679 for(int l=0; l<4; l++){
680 temp_PDF += E[w][j][k][l]*tt1*tt2*tt3*t1rhop[l]*G[l][l];
681 }
682 }
683 }
684 }
685 if(B_Kst0rhop1<0.0){
686 B_Kst0rhop1 = barrier(1,sD,sKst0,srhop,rD2,mDsM);
687 tmp1 = B_Kst0*B_rhop*B_Kst0rhop1*temp_PDF;
688 }
689 }
690 else if(g2[i]==2){
691 calt2(pKstr0,prhop,t2D);
692 for(int w=0; w<4; w++){
693 tt1=t1Kst0[w]*G[w][w];
694 for(int j=0; j<4; j++){
695 temp_PDF += t2D[w][j]*t1rhop[j]*G[j][j]*tt1;
696 }
697 }
698 if (B_Kst0rhop2<0.0){
699 B_Kst0rhop2 = barrier(2,sD,sKst0,srhop,rD2,mDsM);
700 tmp1 = B_Kst0*B_rhop*B_Kst0rhop2*temp_PDF;
701 }
702 }
703 if(isRhop==0){
704 propagatorGS(mass2sq,mass2[i],width2[i],srhop,spionp,spi0,rRes1,proRhop);//rho as mass2
705 isRhop=1.0;
706 }
707 if(g0[i]==1){ propagatorRBWl1(mass1sq,mass1[i],width1[i],sKst0,skaon,spionm,rRes1,pro0);}
708 else if(g0[i]==0){ pro0[0] = 1; pro0[1] = 0;}
709 if(g1[i]==1){ pro1[0] = proRhop[0]; pro1[1] = proRhop[1];}
710 else if(g1[i]==0){ pro1[0] = 1; pro1[1] = 0;}
711 Com_Multi(pro0,pro1,pro);
712 amp_tmp[0] = tmp1*pro[0];
713 amp_tmp[1] = tmp1*pro[1];
714 }
715 if(modetype[i]==2){
716 if (B_Kstp<0.0) B_Kstp = barrier(1,sKstp,skaon,spi0,rRes1,mass1[i]);
717 if (B_rho0<0.0) B_rho0 = barrier(1,srho0,spionp,spionm,rRes1,mass2[i]);
718 if(g2[i]==0){
719 for(int w=0; w<4; w++){
720 temp_PDF += G[w][w]*t1Kstp[w]*t1rho0[w];
721 }
722 tmp1 = B_Kstp*B_rho0*temp_PDF;
723 }
724 else if(g2[i]==1){
725 calt1(pKstrp,prho0,t1D);
726 for(int w=0; w<4; w++){
727 tt1=pD[w]*G[w][w];
728 for(int j=0; j<4; j++){
729 tt2=t1D[j]*G[j][j];
730 for(int k=0; k<4; k++){
731 tt3=t1Kstp[k]*G[k][k];
732 for(int l=0; l<4; l++){
733 temp_PDF += E[w][j][k][l]*tt1*tt2*tt3*t1rho0[l]*G[l][l];
734 }
735 }
736 }
737 }
738 if(B_Kstprho01<0.0){
739 B_Kstprho01 = barrier(1,sD,sKstp,srho0,rD2,mDsM);
740 tmp1 = B_Kstp*B_rho0*B_Kstprho01*temp_PDF;
741 }
742 }
743 else if(g2[i]==2){
744 calt2(pKstrp,prho0,t2D);
745 for(int w=0; w<4; w++){
746 tt1=t1Kstp[w]*G[w][w];
747 for(int j=0; j<4; j++){
748 temp_PDF += t2D[w][j]*t1rho0[j]*G[j][j]*tt1;
749 }
750 }
751 if (B_Kstprho02<0.0){
752 B_Kstprho02 = barrier(2,sD,sKstp,srho0,rD2,mDsM);
753 tmp1 = B_Kstp*B_rho0*B_Kstprho02*temp_PDF;
754 }
755 }
756 if(isRho0==0){
757 propagatorGS(mass2sq,mass2[i],width2[i],srho0,spionp,spionm,rRes1,proRho0);//rho as mass2
758 isRho0=1.0;
759 }
760 if(g0[i]==1){ propagatorRBWl1(mass1sq,mass1[i],width1[i],sKstp,skaon,spi0,rRes1,pro0);}
761 else if(g0[i]==0){ pro0[0] = 1; pro0[1] = 0;}
762 if(g1[i]==1){ pro1[0] = proRho0[0]; pro1[1] = proRho0[1];}
763 else if(g1[i]==0){ pro1[0] = 1; pro1[1] = 0;}
764 Com_Multi(pro0,pro1,pro);
765 amp_tmp[0] = tmp1*pro[0];
766 amp_tmp[1] = tmp1*pro[1];
767 }
768else if(modetype[i]==6){
769 if(B_D_A1<0.0) B_D_A1 = barrier(1,sD,sa1,skaon,rD2,mDsM);
770 if(B_rhop<0.0) B_rhop = barrier(1,srhop,spionp,spi0,rRes1,mass2[i]);
771 calt1(pA1,Kp,t1D);
772 if(g2[i]==0){
773 for(int w=0; w<4; w++){
774 tt1 = t1D[w]*G[w][w];
775 for(int j=0; j<4; j++){
776 temp_PDF += tt1*(G[w][j]-pA1[w]*pA1[j]/sa1)*t1rhop[j]*G[j][j];
777 }
778 }
779 tmp1 = B_rhop*B_D_A1*temp_PDF;
780 }
781 else if(g2[i]==2){
782 for(int w=0; w<4; w++){
783 tt1=t1D[w]*G[w][w];
784 for(int j=0; j<4; j++){
785 temp_PDF += tt1*t2A1[w][j]*t1rhop[j]*G[j][j];
786 }
787 }
788 if(B_D_A2<0.0) B_D_A2 = barrier(2,sa1,srhop,spionm,rRes1,mass1[i]);
789 tmp1 = B_rhop*B_D_A2*B_D_A1*temp_PDF;
790 }
791 if(isRhop==0){
792 propagatorGS(mass2sq,mass2[i],width2[i],srhop,spionp,spi0,rRes1,proRhop);
793 isRhop=1;
794 }
795 if(g0[i]==1) propagatorRBW(mass1sq,mass1[i],width1[i],sa1,srhop,spionm,rRes1,g2[i],pro0);
796 else if(g0[i]==0){pro0[0] = 1; pro0[1] = 0;}
797 if(g1[i]==1){pro1[0] = proRhop[0]; pro1[1] = proRhop[1];}
798 else if(g1[i]==0){pro1[0] = 1; pro1[1] = 0;}
799 Com_Multi(pro0,pro1,pro);
800 amp_tmp[0] = tmp1*pro[0];
801 amp_tmp[1] = tmp1*pro[1];
802 }
803else if(modetype[i]==7){
804 if(B_D_A1<0.0) B_D_A1 = barrier(1,sD,sa1,skaon,rD2,mDsM);
805 if(B_rhom<0.0) B_rhom = barrier(1,srhom,spionm,spi0,rRes1,mass2[i]);
806 calt1(pA1,Kp,t1D);
807 if(g2[i]==0){
808 for(int w=0; w<4; w++){
809 tt1 = t1D[w]*G[w][w];
810 for(int j=0; j<4; j++){
811 temp_PDF += tt1*(G[w][j]-pA1[w]*pA1[j]/sa1)*t1rhom[j]*G[j][j];
812 }
813 }
814 tmp1 = B_rhom*B_D_A1*temp_PDF;
815 }
816 else if(g2[i]==2){
817 for(int w=0; w<4; w++){
818 tt1=t1D[w]*G[w][w];
819 for(int j=0; j<4; j++){
820 temp_PDF += tt1*t2A2[w][j]*t1rhom[j]*G[j][j];
821 }
822 }
823 if(B_D_A2<0.0) B_D_A2 = barrier(2,sa1,srhom,spionp,rRes1,mass1[i]);
824 tmp1 = B_rhom*B_D_A2*B_D_A1*temp_PDF;
825 }
826 if(isRhom==0){
827 propagatorGS(mass2sq,mass2[i],width2[i],srhom,spionm,spi0,rRes1,proRhom);
828 isRhom=1;
829 }
830 if(g0[i]==1) propagatorRBW(mass1sq,mass1[i],width1[i],sa1,srhom,spionp,rRes1,g2[i],pro0);
831 else if(g0[i]==0){pro0[0] = 1; pro0[1] = 0;}
832 if(g1[i]==1){pro1[0] = proRhom[0]; pro1[1] = proRhom[1];}
833 else if(g1[i]==0){pro1[0] = 1; pro1[1] = 0;}
834 Com_Multi(pro0,pro1,pro);
835 amp_tmp[0] = -tmp1*pro[0];
836 amp_tmp[1] = -tmp1*pro[1];
837
838 }
839
840 else if(modetype[i]==10){
841 if(B_rhom<0.0) B_rhom = barrier(1,srhom,spionm,spi0,rRes1,mass2[i]);
842 //printf("B_rhom= %.15f\n",B_rhom);
843 if(B_D_K1rho<0.0) B_D_K1rho = barrier(1,sD,sk11,spionp,rD2,mDsM);
844 //printf("B_D_K1rho= %.15f\n",B_D_K1rho);
845 calt1(pK11,Pip,t1D);
846 if(g2[i]==0){
847 for(int w=0; w<4; w++){
848 tt1=t1D[w]*G[w][w];
849 for(int j=0; j<4; j++){
850 temp_PDF += tt1*(G[w][j]-pK11[w]*pK11[j]/sk11)*t1rhom[j]*G[j][j];
851 }
852 }
853 tmp1 = B_rhom*B_D_K1rho*temp_PDF;
854 }
855 else if(g2[i]==2){
856 for(int w=0; w<4; w++){
857 tt1=t1D[w]*G[w][w];
858 for(int j=0; j<4; j++){
859 temp_PDF += tt1*t2k21[w][j]*t1rhom[j]*G[j][j];
860 }
861 }
862 if(B_K1rho<0.0){
863 B_K1rho = barrier(2,sk11,srhom,skaon,rRes1,mass1[i]);
864 tmp1 = B_rhom*B_K1rho*B_D_K1rho*temp_PDF;
865 }
866 }
867 if(isRhom==0){
868 propagatorGS(mass2sq,mass2[i],width2[i],srhom,spionm,spi0,rRes1,proRhom);//rho also as mass2, but for this process the input parameter should change position (differ to bofore)).
869 isRhom=1.0;
870 }
871 qqq = 0.25*(sk11+srhom-skaon)*(sk11+srhom-skaon)/sk11-srhom;
872 qqq0 = 0.25*(mK1270+srhom-skaon)*(mK1270+srhom-skaon)/mK1270-srhom;
873 if(qqq<2e-4)continue;
874 if(g0[i]==1) propagatorRBW(mass1sq,mass1[i],width1[i],sk11,srhom,skaon,rRes1,g2[i],pro0);//K1270(rho) as mass1 which is different from K1270(K*)(as mass2)
875 else if(g0[i]==0){pro0[0] = 1; pro0[1] = 0;}
876 if(g1[i]==1){pro1[0] = proRhom[0]; pro1[1] = proRhom[1];}
877 else if(g1[i]==0){pro1[0] = 1; pro1[1] = 0;}
878 Com_Multi(pro0,pro1,pro);
879 amp_tmp[0] = tmp1*pro[0];
880 amp_tmp[1] = tmp1*pro[1];
881
882 }
883else if(modetype[i]==11){
884 if(B_D_K1Pi<0.0) B_D_K1Pi = barrier(1,sD,sk11,spionp,rD2,mDsM);
885 if(B_Kst0<0.0) B_Kst0 = barrier(1,sKst0,skaon,spionm,rRes1,mass1[i]);
886 calt1(pK11,Pip,t1D);
887 if(g2[i]==0){
888 for(int w=0; w<4; w++){
889 tt1 = t1D[w]*G[w][w];
890 for(int j=0; j<4; j++){
891 temp_PDF += tt1*(G[w][j]-pK11[w]*pK11[j]/sk11)*t1Kst0[j]*G[j][j];
892 }
893 }
894 tmp1 = B_Kst0*B_D_K1Pi*temp_PDF;
895 }
896 else if(g2[i]==2){
897 for(int w=0; w<4; w++){
898 tt1=t1D[w]*G[w][w];
899 for(int j=0; j<4; j++){
900 temp_PDF += tt1*t2k11[w][j]*t1Kst0[j]*G[j][j];
901 }
902 }
903 if(B_K1rhop2<0.0) B_K1rhop2 = barrier(2,sk11,sKst0,spi0,rRes1,mass2[i]);
904 tmp1 = B_Kst0*B_K1rhop2*B_D_K1Pi*temp_PDF;
905 }
906 if(isKst0==0){
907 propagatorRBWl1(mass1sq,mass1[i],width1[i],sKst0,skaon,spionm,rRes1,proKst0);//Kst0 also as mass1
908 isKst0=1;
909 }
910 if(g0[i]==1){
911 pro0[0] = proKst0[0]; pro0[1] = proKst0[1];
912 }
913 else if(g0[i]==0){pro0[0] = 1; pro0[1] = 0;}
914 if(g1[i]==1){
915 propagatorRBW(mass2sq,mass2[i],width2[i],sk11,sKst0,spi0,rRes1,g2[i],pro1);//K1270(K*) as mass2
916 }
917 else if(g1[i]==0){pro1[0] = 1; pro1[1] = 0;}
918 Com_Multi(pro0,pro1,pro);
919 amp_tmp[0] = tmp1*pro[0];
920 amp_tmp[1] = tmp1*pro[1];
921
922 }
923 else if(modetype[i]==12){
924 if(B_D_K1Pi<0.0) B_D_K1Pi = barrier(1,sD,sk11,spionp,rD2,mDsM);
925 if(B_Kstp<0.0) B_Kstp = barrier(1,sKstp,skaon,spi0,rRes1,mass1[i]);
926 calt1(pK11,Pip,t1D);
927 if(g2[i]==0){
928 for(int w=0; w<4; w++){
929 tt1 = t1D[w]*G[w][w];
930 for(int j=0; j<4; j++){
931 temp_PDF1 += tt1*(G[w][j]-pK11[w]*pK11[j]/sk11)*t1Kstp[j]*G[j][j];
932 }
933 }
934 tmp1 = B_Kstp*B_D_K1Pi*temp_PDF1;
935 }
936 else if(g2[i]==2){
937 for(int w=0; w<4; w++){
938 tt1=t1D[w]*G[w][w];
939 for(int j=0; j<4; j++){
940 temp_PDF1 += tt1*t2k12[w][j]*t1Kstp[j]*G[j][j];
941 }
942 }
943 if(B_K1Kstp2<0.0) B_K1Kstp2 = barrier(2,sk11,sKstp,spionm,rRes1,mass2[i]);
944 tmp1 = B_Kstp*B_K1Kstp2*B_D_K1Pi*temp_PDF1;
945 }
946 if(isKstp==0){
947 propagatorRBWl1(mass1sq,mass1[i],width1[i],sKstp,skaon,spi0,rRes1,proKstp);//Kst0 also as mass1
948 isKstp=1;
949 }
950 if(g0[i]==1){
951 pro0[0] = proKstp[0]; pro0[1] = proKstp[1];
952 }
953 else if(g0[i]==0){pro0[0] = 1; pro0[1] = 0;}
954 if(g1[i]==1){
955 propagatorRBW(mass2sq,mass2[i],width2[i],sk11,sKstp,spionm,rRes1,g2[i],pro1);//K1270(K*) as mass2
956 }
957 else if(g1[i]==0){pro1[0] = 1; pro1[1] = 0;}
958 Com_Multi(pro0,pro1,pro);
959 amp_tmp[0] = -tmp1*pro[0];
960 amp_tmp[1] = -tmp1*pro[1];
961
962 }
963 else if(modetype[i]==14){
964 double B1_omega1=barrier(1,somega,srhop,spionm,rRes1,mass1[i]);
965 double B1_omega2=barrier(1,somega,srho0,spi0,rRes1,mass1[i]);
966 double B1_omega3=barrier(1,somega,srhom,spionp,rRes1,mass1[i]);
967 calt1(pomega,Kp,t1Ds_omega);
968 double B1_1V11=barrier(1,srhop,spionp,spi0,rRes1,mrho);
969 double B1_1V12=barrier(1,srho0,spionp,spionm,rRes1,mrho0);
970 double B1_1V13=barrier(1,srhom,spi0,spionm,rRes1,mrho);
971 double B1_Ds_omega=barrier(1,sD,somega,skaon,rD2,mDsM);
972 for(int w=0; w<4; w++){
973 tt1=pomega[w]*G[w][w];
974 for(int j=0; j<4; j++){
975 tt21=(prhop[j]-Pim[j])*G[j][j];
976 tt22=(prho0[j]-Pi0[j])*G[j][j];
977 tt23=(prhom[j]-Pip[j])*G[j][j];
978 for(int k=0; k<4; k++){
979 tt3=Kp[k]*G[k][k];
980 for(int l=0; l<4; l++){
981 temp_PDF11 += E[w][j][k][l]*tt1*tt21*tt3*(Pip[l]-Pi0[l])*G[l][l];
982 temp_PDF12 += E[w][j][k][l]*tt1*tt22*tt3*(Pip[l]-Pim[l])*G[l][l];
983 temp_PDF13 += E[w][j][k][l]*tt1*tt23*tt3*(Pim[l]-Pi0[l])*G[l][l];
984 }
985 }
986 }
987 }
988 temp_PDF11 = temp_PDF11*B1_Ds_omega*B1_omega1*B1_1V11;
989 temp_PDF12 = temp_PDF12*B1_Ds_omega*B1_omega2*B1_1V12;
990 temp_PDF13 = temp_PDF13*B1_Ds_omega*B1_omega3*B1_1V13;
991 double pro1V11[2],pro1V12[2],pro1V13[2];
992 propagatorGS(mass2sq,mass2[i],width2[i],srhop,spionp,spi0,rRes1,pro1V11);
993 propagatorGS(mass2sq,mass2[i],width2[i],srho0,spionp,spionm,rRes1,pro1V12);
994 propagatorGS(mass2sq,mass2[i],width2[i],srhom,spi0,spionm,rRes1,pro1V13);
995 if(g0[i]==1){
996 propagatorRBW(mass1sq,mass1[i],width1[i],somega,srhop,spionm,rRes1,1,pro0_11);
997 propagatorRBW(mass1sq,mass1[i],width1[i],somega,srho0,spi0,rRes1,1,pro0_12);
998 propagatorRBW(mass1sq,mass1[i],width1[i],somega,srhom,spionp,rRes1,1,pro0_13);
999 }
1000 else if(g0[i]==0){
1001 pro0_11[0] = 1; pro0_11[1] = 0;
1002 pro0_12[0] = 1; pro0_12[1] = 0;
1003 pro0_13[0] = 1; pro0_13[1] = 0;
1004 }
1005 Com_Multi(pro0_11,pro1V11,pro_11);
1006 Com_Multi(pro0_12,pro1V12,pro_12);
1007 Com_Multi(pro0_13,pro1V13,pro_13);
1008 amp_tmp[0] = (temp_PDF11*pro_11[0] + temp_PDF12*pro_12[0] + temp_PDF13*pro_13[0]);
1009 amp_tmp[1] = (temp_PDF11*pro_11[1] + temp_PDF12*pro_12[1] + temp_PDF13*pro_13[1]);
1010
1011 }
1012if(modetype[i]==15){
1013 tmp1 = 1;
1014 if(isPipPi0_S==0){
1015 proPiPi_S[0] = realpipis;
1016 proPiPi_S[1] = imagpipis;
1017 isPipPi0_S=1;
1018 }
1019 if(isKpPim_S==0){
1020 KPiSLASS(sKstp,skaon,spi0,proKPi_S);
1021 isKpPim_S=1;
1022 }
1023 Com_Multi(proKPi_S,proPiPi_S,amp_tmp);
1024
1025 }
1026else if(modetype[i]==16){
1027 if(B_D_A1<0.0) B_D_A1 = barrier(1,sD,sa1,skaon,rD2,mDsM);
1028 if(B_rhop<0.0) B_rhop = barrier(1,srhop,spionp,spi0,rRes1,mass2[i]);
1029 calt1(pA1,Kp,t1D);
1030 if(g2[i]==0){
1031 for(int w=0; w<4; w++){
1032 tt1 = t1D[w]*G[w][w];
1033 for(int j=0; j<4; j++){
1034 temp_PDF += tt1*(G[w][j]-pA1[w]*pA1[j]/sa1)*t1rhop[j]*G[j][j];
1035 }
1036 }
1037 tmp1 = B_rhop*B_D_A1*temp_PDF;
1038 }
1039 else if(g2[i]==2){
1040 for(int w=0; w<4; w++){
1041 tt1=t1D[w]*G[w][w];
1042 for(int j=0; j<4; j++){
1043 temp_PDF += tt1*t2A1[w][j]*t1rhop[j]*G[j][j];
1044 }
1045 }
1046 if(B_D_A2<0.0) B_D_A2 = barrier(2,sa1,srhop,spionm,rRes1,mass1[i]);
1047 tmp1 = B_rhop*B_D_A2*B_D_A1*temp_PDF;
1048 }
1049 if(isRhop==0){
1050 propagatorGS(mass2sq,mass2[i],width2[i],srhop,spionp,spi0,rRes1,proRhop);
1051 isRhop=1;
1052 }
1053 if(g0[i]==1) propagatorRBW(mass1sq,mass1[i],width1[i],sa1,srhop,spionm,rRes1,g2[i],pro0);
1054 else if(g0[i]==0){pro0[0] = 1; pro0[1] = 0;}
1055 if(g1[i]==1){pro1[0] = proRhop[0]; pro1[1] = proRhop[1];}
1056 else if(g1[i]==0){pro1[0] = 1; pro1[1] = 0;}
1057 Com_Multi(pro0,pro1,pro);
1058 amp_tmp[0] = tmp1*pro[0];
1059 amp_tmp[1] = tmp1*pro[1];
1060
1061 }
1062else if(modetype[i]==17){
1063 if(B_D_A1<0.0) B_D_A1 = barrier(1,sD,sa1,skaon,rD2,mDsM);
1064 if(B_rhom<0.0) B_rhom = barrier(1,srhom,spionm,spi0,rRes1,mass2[i]);
1065 calt1(pA1,Kp,t1D);
1066 if(g2[i]==0){
1067 for(int w=0; w<4; w++){
1068 tt1 = t1D[w]*G[w][w];
1069 for(int j=0; j<4; j++){
1070 temp_PDF += tt1*(G[w][j]-pA1[w]*pA1[j]/sa1)*t1rhom[j]*G[j][j];
1071 }
1072 }
1073 tmp1 = B_rhom*B_D_A1*temp_PDF;
1074 }
1075 else if(g2[i]==2){
1076 for(int w=0; w<4; w++){
1077 tt1=t1D[w]*G[w][w];
1078 for(int j=0; j<4; j++){
1079 temp_PDF += tt1*t2A2[w][j]*t1rhom[j]*G[j][j];
1080 }
1081 }
1082 if(B_D_A2<0.0) B_D_A2 = barrier(2,sa1,srhom,spionp,rRes1,mass1[i]);
1083 tmp1 = B_rhom*B_D_A2*B_D_A1*temp_PDF;
1084 }
1085 if(isRhom==0){
1086 propagatorGS(mass2sq,mass2[i],width2[i],srhom,spionm,spi0,rRes1,proRhom);
1087 isRhom=1;
1088 }
1089 if(g0[i]==1) propagatorRBW(mass1sq,mass1[i],width1[i],sa1,srhom,spionp,rRes1,g2[i],pro0);
1090 else if(g0[i]==0){pro0[0] = 1; pro0[1] = 0;}
1091 if(g1[i]==1){pro1[0] = proRhom[0]; pro1[1] = proRhom[1];}
1092 else if(g1[i]==0){pro1[0] = 1; pro1[1] = 0;}
1093 Com_Multi(pro0,pro1,pro);
1094 amp_tmp[0] = -tmp1*pro[0];
1095 amp_tmp[1] = -tmp1*pro[1];
1096 }
1097 Com_Multi(amp_tmp,cof,amp_PDF);
1098 PDF[0] += amp_PDF[0];
1099 PDF[1] += amp_PDF[1];
1100 double value = PDF[0]*PDF[0] + PDF[1]*PDF[1];
1101 if(value <=0) {value = 1e-20;}
1102 Result = value;
1103 q1270 = qqq;
1104}
1105}
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
s
XmlRpcServer s
Definition HelloServer.cpp:11
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
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
EvtDsToKpPipPimPi0::clone
EvtDecayBase * clone()
Definition EvtDsToKpPipPimPi0.cc:43
EvtDsToKpPipPimPi0::decay
void decay(EvtParticle *p)
Definition EvtDsToKpPipPimPi0.cc:161
EvtDsToKpPipPimPi0::getName
void getName(std::string &name)
Definition EvtDsToKpPipPimPi0.cc:39
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
m2
double double * m2
Definition qcdloop1.h:75
b
const double b
Definition slope.cxx:9