BOSS 7.1.0
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtDToKSpipipi.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: EvtDToKSpipipi.cc
11// the necessary file: EvtDToKSpipipi.hh
12//
13// Description: D+ -> KS0 pi+ pi+ pi-,
14// PHYSICAL REVIEW D 100, 072008 (2019)
15//
16// Modification history:
17//
18// Liaoyuan Dong Jan.15, 2020 Module created
19//
20//------------------------------------------------------------------------
21#include "EvtGenBase/EvtPatches.hh"
22#include "EvtGenBase/EvtParticle.hh"
23#include "EvtGenBase/EvtGenKine.hh"
24#include "EvtGenBase/EvtPDL.hh"
25#include "EvtGenBase/EvtReport.hh"
26#include "EvtGenModels/EvtDToKSpipipi.hh"
27#include "EvtGenBase/EvtComplex.hh"
28#include "EvtGenBase/EvtDecayTable.hh"
29#include <stdlib.h>
30
31EvtDToKSpipipi::~EvtDToKSpipipi() {}
32
33void EvtDToKSpipipi::getName(std::string& model_name){
34 model_name="DToKSpipipi";
35}
36
37EvtDecayBase* EvtDToKSpipipi::clone(){
38 return new EvtDToKSpipipi;
39}
40
41void EvtDToKSpipipi::init(){
42 checkNArg(0);
43 checkNDaug(4);
44 checkSpinParent(EvtSpinType::SCALAR);
45/*
46 checkSpinDaughter(0,EvtSpinType::SCALAR);
47 checkSpinDaughter(1,EvtSpinType::SCALAR);
48 checkSpinDaughter(3,EvtSpinType::SCALAR);
49 checkSpinDaughter(4,EvtSpinType::SCALAR);
50*/
51 //std::cout << "Initializing EvtDToKSpipipi" << std::endl;
52
53 mK1270 = 1.272; mK1400 = 1.403;
54 GK1270 = 0.09; GK1400 = 0.174;
55 mKstr = 0.89166; mrho = 0.77549;
56 GKstr = 0.0462; Grho = 0.1491;
57 msigma = 0.472;
58 Gsigma = 0.542;
59 phi_omega = -0.02;
60 mK1650 = 1.65;
61 GK1650 = 0.158;
62 rho[0] = 1.0;
63 phi[0] = 0.0;
64
65 ma1 = 1.22;
66 Ga1 = 0.4282;
67 mK1460 = 1.4152;
68 GK1460 = 0.2485;
69 rho_omega = 0.00294;
70
71 phi[1] = -1.55;
72 rho[1] = 0.5843;
73
74 phi[2] = -1.8223;
75 rho[2] = 2.0974;
76
77 phi[3] = -2.6751;
78 rho[3] = 0.46642;
79
80 phi[4] = -2.2429;
81 rho[4] = 0.33334;
82
83 phi[5] = -0.55888;
84 rho[5] = 0.15549;
85
86 phi[6] = -1.8778;
87 rho[6] = 0.94452;
88
89 phi[7] = 2.7724;
90 rho[7] = 0.99411;
91
92 phi[8] = -2.6461;
93 rho[8] = 0.21231;
94
95 phi[9] = -0.95137;
96 rho[9] = 0.29895;
97
98 phi[10] = -3.0843;
99 rho[10] = 3.6361;
100
101 phi[11] = 2.0954;
102 rho[11] = 0.96472;
103
104 phi[12] = -2.4965;
105 rho[12] = 0.48470;
106
107 //for (int i=0; i<13; i++) {
108 // cout << i << " rho,phi = " << rho[i] << ", "<< phi[i] << endl;
109 //}
110
111 mD = 1.86486;
112 rD = 5;
113 metap = 0.95778;
114 mkstr = 0.89594;
115 mk0 = 0.497614;
116 mass_Kaon = 0.49368;
117 mass_Pion = 0.13957;
118 math_pi = 3.1415926;
119
120 pi = 3.1415926;
121 mpi = 0.13957;
122 g1 = 0.5468;
123 g2 = 0.23;
124
125 int GG[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
126 int EE[4][4][4][4] =
127 { { {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
128 {{0,0,0,0}, {0,0,0,0}, {0,0,0,1}, {0,0,-1,0}},
129 {{0,0,0,0}, {0,0,0,-1}, {0,0,0,0}, {0,1,0,0} },
130 {{0,0,0,0}, {0,0,1,0}, {0,-1,0,0}, {0,0,0,0} } },
131 { {{0,0,0,0}, {0,0,0,0}, {0,0,0,-1}, {0,0,1,0} },
132 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
133 {{0,0,0,1}, {0,0,0,0}, {0,0,0,0}, {-1,0,0,0}},
134 {{0,0,-1,0}, {0,0,0,0}, {1,0,0,0}, {0,0,0,0} } },
135 { {{0,0,0,0}, {0,0,0,1}, {0,0,0,0}, {0,-1,0,0}},
136 {{0,0,0,-1}, {0,0,0,0}, {0,0,0,0}, {1,0,0,0} },
137 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
138 {{0,1,0,0}, {-1,0,0,0}, {0,0,0,0}, {0,0,0,0} } },
139 { {{0,0,0,0}, {0,0,-1,0}, {0,1,0,0}, {0,0,0,0} },
140 {{0,0,1,0}, {0,0,0,0}, {-1,0,0,0}, {0,0,0,0} },
141 {{0,-1,0,0}, {1,0,0,0}, {0,0,0,0}, {0,0,0,0} },
142 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} } } };
143 for (int i=0; i<4; i++) {
144 for (int j=0; j<4; j++) {
145 G[i][j] = GG[i][j];
146 for (int k=0; k<4; k++) {
147 for (int l=0; l<4; l++) {
148 E[i][j][k][l] = EE[i][j][k][l];
149 }
150 }
151 }
152 }
153}
154
155void EvtDToKSpipipi::initProbMax() {
156 setProbMax(3700.0);
157}
158
159void EvtDToKSpipipi::decay( EvtParticle *p){
160/*
161 double maxprob = 0.0;
162 for(int ir=0;ir<=60000000;ir++){
163 p->initializePhaseSpace(getNDaug(),getDaugs());
164 EvtVector4R Ks0 = p->getDaug(0)->getP4();
165 EvtVector4R pi1 = p->getDaug(1)->getP4();
166 EvtVector4R pi2 = p->getDaug(2)->getP4();
167 EvtVector4R pi3 = p->getDaug(3)->getP4();
168 double Ks[4],Pip1[4],Pip2[4],Pim[4];
169 Ks[0] = Ks0.get(0); Pip1[0] = pi1.get(0); Pip2[0] = pi2.get(0); Pim[0] = pi3.get(0);
170 Ks[1] = Ks0.get(1); Pip1[1] = pi1.get(1); Pip2[1] = pi2.get(1); Pim[1] = pi3.get(1);
171 Ks[2] = Ks0.get(2); Pip1[2] = pi1.get(2); Pip2[2] = pi2.get(2); Pim[2] = pi3.get(2);
172 Ks[3] = Ks0.get(3); Pip1[3] = pi1.get(3); Pip2[3] = pi2.get(3); Pim[3] = pi3.get(3);
173 double Prob = calPDF(Ks, Pip1, Pip2, Pim);
174 if(Prob>maxprob) {
175 maxprob=Prob;
176 std::cout << "Max PDF = " << ir << " prob= " << Prob << std::endl;
177 }
178 }
179 std::cout << "Max!!!!!!!!!!! " << maxprob<< std::endl;
180*/
181 p->initializePhaseSpace(getNDaug(),getDaugs());
182 EvtVector4R Ks0 = p->getDaug(0)->getP4();
183 EvtVector4R pi1 = p->getDaug(1)->getP4();
184 EvtVector4R pi2 = p->getDaug(2)->getP4();
185 EvtVector4R pi3 = p->getDaug(3)->getP4();
186
187 double Ks[4],Pip1[4],Pip2[4],Pim[4];
188 Ks[0] = Ks0.get(0); Pip1[0] = pi1.get(0); Pip2[0] = pi2.get(0); Pim[0] = pi3.get(0);
189 Ks[1] = Ks0.get(1); Pip1[1] = pi1.get(1); Pip2[1] = pi2.get(1); Pim[1] = pi3.get(1);
190 Ks[2] = Ks0.get(2); Pip1[2] = pi1.get(2); Pip2[2] = pi2.get(2); Pim[2] = pi3.get(2);
191 Ks[3] = Ks0.get(3); Pip1[3] = pi1.get(3); Pip2[3] = pi2.get(3); Pim[3] = pi3.get(3);
192 double prob = calPDF(Ks, Pip1, Pip2, Pim);
193 setProb(prob);
194 return;
195}
196
197double EvtDToKSpipipi::calPDF(double Ks[], double Pip1[], double Pip2[], double Pim[]) {
198 EvtComplex PDF[100];
199 double P14[4], P24[4], P34[4], P124[4], P134[4], P234[4];
200 for(int i=0; i<4; i++){
201 P14[i] = Ks[i] + Pim[i];
202 P24[i] = Pip1[i] + Pim[i];
203 P34[i] = Pip2[i] + Pim[i];
204 P124[i] = P14[i] + Pip1[i];
205 P134[i] = P14[i] + Pip2[i];
206 P234[i] = P24[i] + Pip2[i];
207 }
208 //----------D->a1Ks--------------
209 //----------a1->rhoPi------------
210 PDF[0] = D2AP_A2VP(Ks,Pip2,Pip1,Pim,0)*getprop(P24,Pip2,ma1,Ga1,1,0)*
211 getprop(Pip1,Pim,mrho,Grho,2,1) +
212 D2AP_A2VP(Ks,Pip1,Pip2,Pim,0)*getprop(P34,Pip1,ma1,Ga1,1,0)*
213 getprop(Pip2,Pim,mrho,Grho,2,1);
214 PDF[1] = D2AP_A2VP(Ks,Pip2,Pip1,Pim,2)*getprop(P24,Pip2,ma1,Ga1,1,2)*
215 getprop(Pip1,Pim,mrho,Grho,2,1) +
216 D2AP_A2VP(Ks,Pip1,Pip2,Pim,2)*getprop(P34,Pip1,ma1,Ga1,1,2)*
217 getprop(Pip2,Pim,mrho,Grho,2,1);
218 //----------a1->sigma pi---------
219 PDF[2] = D2AP_A2SP(Ks,Pip2,Pip1,Pim)*getprop(P24,Pip2,ma1,Ga1,1,1)*
220 getprop(Pip1,Pim,msigma,Gsigma,4,0) +
221 D2AP_A2SP(Ks,Pip1,Pip2,Pim)*getprop(P34,Pip1,ma1,Ga1,1,1)*
222 getprop(Pip2,Pim,msigma,Gsigma,4,0);
223 //----------D->a1K finish-----
224 //---------D->K1(1400) pi-----
225 //K1400[S]->K* pi
226 PDF[3] = D2AP_A2VP(Pip2,Pip1,Ks,Pim,0)*getprop(P14,Pip1,mK1400,GK1400,1,0)*
227 getprop(Ks,Pim,mKstr,GKstr,1,1) +
228 D2AP_A2VP(Pip1,Pip2,Ks,Pim,0)*getprop(P14,Pip2,mK1400,GK1400,1,0)*
229 getprop(Ks,Pim,mKstr,GKstr,1,1);
230 //K1400[D]->K* pi
231 PDF[4] = D2AP_A2VP(Pip2,Pip1,Ks,Pim,2)*getprop(P14,Pip1,mK1400,GK1400,1,2)*
232 getprop(Ks,Pim,mKstr,GKstr,1,1) +
233 D2AP_A2VP(Pip1,Pip2,Ks,Pim,2)*getprop(P14,Pip2,mK1400,GK1400,1,2)*
234 getprop(Ks,Pim,mKstr,GKstr,1,1);
235 //-----------------------------
236 //-------K1270[S]->Ksrho-------
237 PDF[5] = D2AP_A2VP(Pip2,Ks,Pip1,Pim,0)*getprop(P24,Ks,mK1270,GK1270,0,0)*
238 getprop(Pip1,Pim,mrho,Grho,2,1) +
239 D2AP_A2VP(Pip1,Ks,Pip2,Pim,0)*getprop(P34,Ks,mK1270,GK1270,0,0)*
240 getprop(Pip2,Pim,mrho,Grho,2,1);
241 //-------D->rhoKAD------------
242 PDF[6] = D2AP_A2VP(Pip2,Ks,Pip1,Pim,0)*getprop(Pip1,Pim,mrho,Grho,2,1) +
243 D2AP_A2VP(Pip1,Ks,Pip2,Pim,0)*getprop(Pip2,Pim,mrho,Grho,2,1);
244 PDF[7] = D2AP_A2VP(Pip2,Ks,Pip1,Pim,2)*getprop(Pip1,Pim,mrho,Grho,2,1) +
245 D2AP_A2VP(Pip1,Ks,Pip2,Pim,2)*getprop(Pip2,Pim,mrho,Grho,2,1);
246 //-------D->K1460, K1460->Ks rho---------
247 PDF[8] = D2PP_P2VP(Pip2,Ks,Pip1,Pim)*getprop(P24,Ks,mK1460,GK1460,1,1)*
248 getprop(Pip1,Pim,mrho,Grho,2,1) +
249 D2PP_P2VP(Pip1,Ks,Pip2,Pim)*getprop(P34,Ks,mK1460,GK1460,1,1)*
250 getprop(Pip2,Pim,mrho,Grho,2,1);
251 //--------K*PiA (K1650)---------------------
252 PDF[9] = D2AP_A2VP(Pip2,Pip1,Ks,Pim,0)*getprop(P14,Pip1,mK1650,GK1650,1,0)*
253 getprop(Ks,Pim,mKstr,GKstr,1,1) +
254 D2AP_A2VP(Pip1,Pip2,Ks,Pim,0)*getprop(P14,Pip2,mK1650,GK1650,1,0)*
255 getprop(Ks,Pim,mKstr,GKstr,1,1);
256 //-------KsPiPiSPiA-----------------
257 PDF[10] = D2AP_A2SP(Pip2,Ks,Pip1,Pim) + D2AP_A2SP(Pip1,Ks,Pip2,Pim);
258 //-------KPiS wave------------------
259 EvtComplex corr(2,0);
260 PDF[11] = corr*PHSP(Ks,Pim);
261 // cout<<PHSP(Ks,Pim)<<endl;
262 //-------D->K1460pi, K1460->K*-pi+--------
263 PDF[12] = D2PP_P2VP(Pip2,Pip1,Ks,Pim)*getprop(P14,Pip1,mK1460,GK1460,1,1)*
264 getprop(Ks,Pim,mKstr,GKstr,1,1) +
265 D2PP_P2VP(Pip1,Pip2,Ks,Pim)*getprop(P14,Pip2,mK1460,GK1460,1,1)*
266 getprop(Ks,Pim,mKstr,GKstr,1,1);
267 //-------------------------------------------
268 EvtComplex cof;
269 EvtComplex pdf, module;
270 pdf = EvtComplex (0,0);
271 for(int i=0; i<13; i++){
272 cof = EvtComplex (rho[i]*cos(phi[i]),rho[i]*sin(phi[i]));
273 pdf = pdf + cof*PDF[i];
274 }
275 module = conj(pdf)*pdf;
276 double value;
277 value = real(module);
278 return (value <= 0) ? 1e-20 : value;
279}
280
281EvtComplex EvtDToKSpipipi::KPiSFormfactor(double sa, double sb, double sc, double r)
282{
283 double m1430 = 1.463;
284 double sa0 = m1430*m1430;
285 double w1430 = 0.233;
286 double q0 = (sa0+sb-sc)*(sa0+sb-sc)/(4*sa0)-sb;
287 if(q0<0) q0 = 1e-16;
288 double qs = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
289 double q = sqrt(qs);
290 double width = w1430*q*m1430/sqrt(sa*q0);
291 double temp_R = atan(m1430*width/(sa0-sa));
292 if(temp_R<0) temp_R += math_pi;
293 double deltaR = -5.31 + temp_R;
294 double temp_F = atan(2*1.07*q/(2+(-1.8)*1.07*qs));
295 if(temp_F<0) temp_F += math_pi;
296 double deltaF = 2.33 + temp_F;
297 EvtComplex cR(cos(deltaR), sin(deltaR));
298 EvtComplex cF(cos(deltaF), sin(deltaF));
299 EvtComplex amp = 0.8*sin(deltaF)*cF + sin(deltaR)*cR*cF*cF;
300 return amp;
301}
302EvtComplex EvtDToKSpipipi::D2AP_A2VP(double P1[], double P2[], double P3[], double P4[], int L)
303{
304 double temp_PDF = 0;
305 EvtComplex amp_PDF(0,0);
306 EvtComplex pro[2];
307 double t1V[4], t1D[4], t2A[4][4];
308 double sa[3], sb[3], sc[3], B[3];
309 double pV[4],pA[4],pD[4];
310 for(int i=0; i!=4; i++){
311 pV[i] = P3[i] + P4[i];
312 pA[i] = pV[i] + P2[i];
313 pD[i] = pA[i] + P1[i];
314 }
315 sa[0] = dot(pV,pV);
316 sb[0] = dot(P3,P3);
317 sc[0] = dot(P4,P4);
318 sa[1] = dot(pA,pA);
319 sb[1] = sa[0];
320 sc[1] = dot(P2,P2);
321 sa[2] = dot(pD,pD);
322 sb[2] = sa[1];
323 sc[2] = dot(P1,P1);
324 B[0] = barrier(1,sa[0],sb[0],sc[0],3);
325 B[2] = barrier(1,sa[2],sb[2],sc[2],rD);
326 calt1(P3,P4,t1V);
327 calt1(pA,P1,t1D);
328 if(L == 0){
329 for(int i=0; i!=4; i++){
330 for(int j=0; j!=4; j++){
331 temp_PDF += t1D[i]*(G[i][j] - pA[i]*pA[j]/sa[1])*t1V[j]*(G[i][i])*(G[j][j]);
332 }
333 }
334 B[1] = 1;
335 }
336 if(L == 2){
337 calt2(pV,P2,t2A);
338 for(int i=0; i!=4; i++){
339 for(int j=0; j!=4; j++){
340 temp_PDF += t1D[i]*t2A[i][j]*t1V[j]*(G[i][i])*(G[j][j]);
341 }
342 }
343 B[1] = barrier(2,sa[1],sb[1],sc[1],3);
344 }
345 amp_PDF = temp_PDF*B[0]*B[1]*B[2];
346 return amp_PDF;
347}
348EvtComplex EvtDToKSpipipi::D2AP_A2SP(double P1[], double P2[], double P3[], double P4[])
349{
350 double temp_PDF = 0;
351 EvtComplex amp_PDF(0,0);
352 EvtComplex pro;
353 double sa[3], sb[3], sc[3], B[3];
354 double t1D[4], t1A[4];
355 double pS[4], pA[4], pD[4];
356 for(int i=0; i!=4; i++){
357 pS[i] = P3[i] + P4[i];
358 pA[i] = pS[i] + P2[i];
359 pD[i] = pA[i] + P1[i];
360 }
361 sa[0] = dot(pS,pS);
362 sb[0] = dot(P3,P3);
363 sc[0] = dot(P4,P4);
364 sa[1] = dot(pA,pA);
365 sb[1] = sa[0];
366 sc[1] = dot(P2,P2);
367 sa[2] = dot(pD,pD);
368 sb[2] = sa[1];
369 sc[2] = dot(P1,P1);
370 B[1] = barrier(1,sa[1],sb[1],sc[1],3);
371 B[2] = barrier(1,sa[2],sb[2],sc[2],rD);
372 calt1(pA,P1,t1D);
373 calt1(pS,P2,t1A);
374 for(int i=0; i!=4; i++){
375 temp_PDF += t1D[i]*t1A[i]*(G[i][i]);
376 }
377 amp_PDF = temp_PDF*B[1]*B[2];
378 return amp_PDF;
379}
380EvtComplex EvtDToKSpipipi::D2PP_P2VP(double P1[], double P2[], double P3[], double P4[])
381{
382 double temp_PDF = 0;
383 EvtComplex amp_PDF(0,0);
384 EvtComplex pro;
385 double sa[3], sb[3], sc[3], B[3];
386 double t1V[4];
387 double pV[4], pP[4], pD[4];
388 for(int i=0; i!=4; i++){
389 pV[i] = P3[i] + P4[i];
390 pP[i] = pV[i] + P2[i];
391 pD[i] = pP[i] + P1[i];
392 }
393 sa[0] = dot(pV,pV);
394 sb[0] = dot(P3,P3);
395 sc[0] = dot(P4,P4);
396 sa[1] = dot(pP,pP);
397 sb[1] = sa[0];
398 sc[1] = dot(P2,P2);
399 sa[2] = dot(pD,pD);
400 sb[2] = sa[1];
401 sc[2] = dot(P1,P1);
402 B[0] = barrier(1,sa[0],sb[0],sc[0],3);
403 B[1] = barrier(1,sa[1],sb[1],sc[1],3);
404 calt1(P3,P4,t1V);
405 for(int i=0; i!=4; i++){
406 temp_PDF += P2[i]*t1V[i]*(G[i][i]);
407 }
408 amp_PDF = temp_PDF*B[0]*B[1];
409 return amp_PDF;
410}
411EvtComplex EvtDToKSpipipi::D2VP_V2VP(double P1[], double P2[], double P3[], double P4[])
412{
413 double temp_PDF = 0;
414 EvtComplex amp_PDF(0,0);
415 EvtComplex pro;
416 double sa[3], sb[3], sc[3], B[3];
417 double pV1[4], pV2[4], qV1[4], qV2[4], pD[4];
418 for(int i=0; i!=4; i++){
419 pV2[i] = P3[i] + P4[i];
420 qV2[i] = P3[i] - P4[i];
421 pV1[i] = pV2[i] + P2[i];
422 qV1[i] = pV2[i] - P2[i];
423 pD[i] = pV1[i] + P1[i];
424 }
425 for(int i=0; i!=4; i++){
426 for(int j=0; j!=4; j++){
427 for(int k=0; k!=4; k++){
428 for(int l=0; l!=4; l++){
429 temp_PDF += E[i][j][k][l]*pV1[i]*qV1[j]*P1[k]*qV2[l]*(G[i][i])*(G[j][j])*(G[k][k])*(G[l][l]);
430 }
431 }
432 }
433 }
434 sa[0] = dot(pV2,pV2);
435 sb[0] = dot(P3,P3);
436 sc[0] = dot(P4,P4);
437 sa[1] = dot(pV1,pV1);
438 sb[1] = sa[0];
439 sc[1] = dot(P2,P2);
440 sa[2] = dot(pD,pD);
441 sb[2] = sa[1];
442 sc[2] = dot(P1,P1);
443 B[0] = barrier(1,sa[0],sb[0],sc[0],3.0);
444 B[1] = barrier(1,sa[1],sb[1],sc[1],3.0);
445 B[2] = barrier(1,sa[2],sb[2],sc[2],rD);
446 amp_PDF = temp_PDF*B[0]*B[1]*B[2];
447 return amp_PDF;
448}
449EvtComplex EvtDToKSpipipi::PHSP(double Km[], double Pip[])
450{
451 EvtComplex amp_PDF(0,0);
452 double sa,sb,sc;
453 double KPi[4];
454 for(int i=0; i!=4; i++){
455 KPi[i] = Km[i] + Pip[i];
456 }
457 sa = dot(KPi,KPi);
458 sb = dot(Km,Km);
459 sc = dot(Pip,Pip);
460 amp_PDF = KPiSFormfactor(sa,sb,sc,3.0);
461 return amp_PDF;
462}
463EvtComplex EvtDToKSpipipi::getprop(double daug1[], double daug2[], double mass, double width, int flag, int L){
464 //flag = 0, RBW with constant width
465 //flag = 1, RBW with width depends on momentum
466 //flag = 2, GS convolute with RBW_omega
467 //flag = 3, flatte for f_0(980)/a_0(980)
468 //flag = 4, Bugg formula for sigma;
469 //effect radii == 3.0 GeV^-1
470 EvtComplex prop1, prop2, prop;
471 EvtComplex one(1,0);
472 double pR[4];
473 for(int i=0; i<4; i++){
474 pR[i] = daug1[i] + daug2[i];
475 }
476 double sa, sb, sc;
477 sa = dot(pR,pR);
478 sb = dot(daug1,daug1);
479 sc = dot(daug2,daug2);
480 if(flag == 0) return propogator(mass,width,sa);
481 if(flag == 1) return propagatorRBW(mass,width,sa,sb,sc,3.0,L);
482 if(flag == 2){
483 prop1 = propagatorGS(mass,width,sa,sb,sc,3.0,L);
484 prop2 = propagatorRBW(0.783,0.008,sa,sb,sc,3.0,L);
485 EvtComplex coef_omega(rho_omega*cos(phi_omega),rho_omega*sin(phi_omega));
486 prop = prop1*(one + 0.783*0.783*coef_omega*prop2);
487 return prop;
488 }
489 if(flag == 3){
490 //Not need for D+ -> Ks 3pi
491 }
492 if(flag == 4){
493 EvtComplex ci(0,1);
494 double f = 0.5843+1.6663*sa;
495 double M = 0.9264;
496 double mpi2 = mass_Pion*mass_Pion;
497 double mass2 = M*M;
498 double g1 = f*(sa-mpi2/2)/(mass2-mpi2/2)*exp((mass2-sa)/1.082);
499 EvtComplex rho1s = rhoab(sa,sb,sc);
500 EvtComplex rho1M = rhoab(mass2,sb,sc);
501 EvtComplex rho2s = rho4Pi(sa);
502 EvtComplex rho2M = rho4Pi(mass2);
503 prop = 1.0/(M*M-sa-ci*M*(g1*rho1s/rho1M+0.0024*rho2s/rho2M));
504 return prop;
505 }
506 cout<<"Please set a right flag! "<<flag<<endl;
507 return one;
508}
509EvtComplex EvtDToKSpipipi::rhoab(double sa, double sb, double sc){
510 EvtComplex one(1,0);
511 double q = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
512 EvtComplex rho;
513 EvtComplex ci(0,1);
514 if(q>0) rho = one*sqrt(q/sa);
515 if(q<0) rho = ci*sqrt(-q/sa);
516 rho = 2.0*rho;
517 return rho;
518}
519EvtComplex EvtDToKSpipipi::rho4Pi(double sa)
520{
521 double mpi = 0.13957;
522 EvtComplex one(1,0);
523 EvtComplex rho(0,0);
524 EvtComplex ci(0,1);
525 double temp = 1-16*mpi*mpi/sa;
526 if(temp > 0) rho = one*sqrt(temp)/(1+exp(9.8-3.5*sa));
527 if(temp < 0) rho = ci*sqrt(-temp)/(1+exp(9.8-3.5*sa));
528 return rho;
529}
530
531EvtComplex EvtDToKSpipipi::propogator(double mass, double width, double sx)const
532{
533 EvtComplex ci(0,1);
534 EvtComplex prop=1.0/(mass*mass-sx-ci*mass*width);
535 return prop;
536}
537double EvtDToKSpipipi::wid(double mass, double sa, double sb, double sc, double r, int l)const
538{
539 double widm(0.), q(0.), q0(0.);
540 double sa0 = mass*mass;
541 double m = sqrt(sa);
542 q = Qabcs(sa,sb,sc);
543 q0 = Qabcs(sa0,sb,sc);
544 double z,z0;
545 z = q*r*r;
546 z0 = q0*r*r;
547 double t = q/q0;
548 double F(0.);
549 if(l == 0) F = 1;
550 if(l == 1) F = sqrt((1+z0)/(1+z));
551 if(l == 2) F = sqrt((9+3*z0+z0*z0)/(9+3*z+z*z));
552 widm = pow(t,l+0.5)*mass/m*F*F;
553 return widm;
554}
555double EvtDToKSpipipi::h(double m, double q)const
556{
557 double h(0.);
558 h = 2/pi*q/m*log((m+2*q)/(2*mpi));
559 return h;
560}
561double EvtDToKSpipipi::dh(double mass, double q0)const
562{
563 double dh = h(mass,q0)*(1.0/(8*q0*q0)-1.0/(2*mass*mass))+1.0/(2*pi*mass*mass);
564 return dh;
565}
566double EvtDToKSpipipi::f(double mass, double sx, double q0, double q)const
567{
568 double m = sqrt(sx);
569 double f = mass*mass/(pow(q0,3))*(q*q*(h(m,q)-h(mass,q0))+(mass*mass-sx)*q0*q0*dh(mass,q0));
570 return f;
571}
572double EvtDToKSpipipi::d(double mass, double q0)const
573{
574 double d = 3.0/pi*mpi*mpi/(q0*q0)*log((mass+2*q0)/(2*mpi))+mass/(2*pi*q0) - (mpi*mpi*mass)/(pi*pow(q0,3));
575 return d;
576}
577EvtComplex EvtDToKSpipipi::propagatorRBW(double mass, double width, double sa, double sb, double sc, double r, int l)const
578{
579 EvtComplex ci(0,1);
580 EvtComplex prop=1.0/(mass*mass-sa-ci*mass*width*wid(mass,sa,sb,sc,r,l));
581 return prop;
582}
583EvtComplex EvtDToKSpipipi::propagatorGS(double mass, double width, double sa, double sb, double sc, double r, int l)const
584{
585 EvtComplex ci(0,1);
586 double q = Qabcs(sa,sb,sc);
587 double sa0 = mass*mass;
588 double q0 = Qabcs(sa0,sb,sc);
589 q = sqrt(q);
590 q0 = sqrt(q0);
591 EvtComplex prop = (1+d(mass,q0)*width/mass)/(mass*mass-sa+width*f(mass,sa,q0,q)-ci*mass*width*wid(mass,sa,sb,sc,r,l));
592 return prop;
593}
594double EvtDToKSpipipi::Flatte_rhoab(double sa, double sb, double sc)const
595{
596 double q = Qabcs(sa,sb,sc);
597 double rho = sqrt(q/sa);
598 return rho;
599}
600EvtComplex EvtDToKSpipipi::propagatorFlatte(double mass, double width, double sx, double *sb, double *sc)const
601{
602 EvtComplex ci(0,1);
603 double rho1 = Flatte_rhoab(sx,sb[0],sc[0]);
604 double rho2 = Flatte_rhoab(sx,sb[1],sc[1]);
605 EvtComplex prop = 1.0/(mass*mass-sx-ci*(g1*g1*rho1+g2*g2*rho2));
606 return prop;
607}
608
609double EvtDToKSpipipi::dot(double *a1, double *a2)const
610{
611 double dot = 0;
612 for(int i=0; i!=4; i++){
613 dot += a1[i]*a2[i]*G[i][i];
614 }
615 return dot;
616}
617double EvtDToKSpipipi::Qabcs(double sa, double sb, double sc)const
618{
619 double Qabcs = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
620 if(Qabcs < 0) Qabcs = 1e-16;
621 return Qabcs;
622}
623double EvtDToKSpipipi::barrier(double l, double sa, double sb, double sc, double r)const
624{
625 double q = Qabcs(sa,sb,sc);
626 q = sqrt(q);
627 double z = q*r;
628 z = z*z;
629 double F = 1;
630 if(l > 2) F = 0;
631 if(l == 0)
632 F = 1;
633 if(l == 1){
634 F = sqrt((2*z)/(1+z));
635 }
636 if(l == 2){
637 F = sqrt((13*z*z)/(9+3*z+z*z));
638 }
639 return F;
640}
641void EvtDToKSpipipi::calt1(double daug1[], double daug2[], double t1[]) const
642{
643 double p, pq;
644 double pa[4], qa[4];
645 for(int i=0; i!=4; i++){
646 pa[i] = daug1[i] + daug2[i];
647 qa[i] = daug1[i] - daug2[i];
648 }
649 p = dot(pa,pa);
650 pq = dot(pa,qa);
651 for(int i=0; i!=4; i++){
652 t1[i] = qa[i] - pq/p*pa[i];
653 }
654}
655void EvtDToKSpipipi::calt2(double daug1[], double daug2[], double t2[][4]) const
656{
657 double p,r;
658 double pa[4], t1[4];
659 calt1(daug1,daug2,t1);
660 r = dot(t1,t1);
661 for(int i=0; i!=4; i++){
662 pa[i] = daug1[i] + daug2[i];
663 }
664 p = dot(pa,pa);
665 for(int i=0; i!=4; i++){
666 for(int j=0; j!=4; j++){
667 t2[i][j] = t1[i]*t1[j] - 1.0/3*r*(G[i][j]-pa[i]*pa[j]/p);
668 }
669 }
670}
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
conj
Evt3Rank3C conj(const Evt3Rank3C &t2)
Definition: Evt3Rank3C.cc:175
exp
EvtComplex exp(const EvtComplex &c)
Definition: EvtComplex.hh:252
mpi
const double mpi
Definition: Gam4pikp.cxx:47
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
t
TTree * t
Definition: binning.cxx:23
EvtDToKSpipipi::getName
void getName(std::string &name)
Definition: EvtDToKSpipipi.cc:33
EvtDToKSpipipi::decay
void decay(EvtParticle *p)
Definition: EvtDToKSpipipi.cc:159
EvtDToKSpipipi::clone
EvtDecayBase * clone()
Definition: EvtDToKSpipipi.cc:37
EvtDToKSpipipi::~EvtDToKSpipipi
virtual ~EvtDToKSpipipi()
Definition: EvtDToKSpipipi.cc:31
EvtDecayBase::checkSpinParent
void checkSpinParent(EvtSpinType::spintype sp)
Definition: EvtDecayBase.cc:520
EvtDecayBase::setProbMax
void setProbMax(double prbmx)
Definition: EvtDecayBase.cc:297
EvtDecayBase::checkNDaug
void checkNDaug(int d1, int d2=-1)
Definition: EvtDecayBase.cc:504
EvtDecayBase::getDaugs
EvtId * getDaugs()
Definition: EvtDecayBase.hh:65
EvtDecayBase::checkNArg
void checkNArg(int a1, int a2=-1, int a3=-1, int a4=-1)
Definition: EvtDecayBase.cc:482
EvtDecayProb::setProb
void setProb(double prob)
Definition: EvtDecayProb.hh:34
EvtParticle::getP4
const EvtVector4R & getP4() const
Definition: EvtParticle.cc:120
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition: EvtParticle.cc:84
EvtParticle::initializePhaseSpace
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
Definition: EvtParticle.cc:1070
EvtVector4R::get
double get(int i) const
Definition: EvtVector4R.hh:179
rb::mpi2
const double mpi2
Definition: TConstant.h:28
one
double precision pisqo6 one
Definition: qlconstants.h:4