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