BOSS 7.0.9
BESIII Offline Software System
Loading...
Searching...
No Matches
BesEvtGen-00-04-08/src/EvtGen/EvtGenModels/EvtD0ToKpipipi.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: EvtD0ToKpipipi.cc
11// the necessary file: EvtD0ToKpipipi.hh
12//
13// Description: D0 -> K- pi+ pi+ pi-,
14// see PHYSICAL REVIEW D 95, 072010 (2017)
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/EvtD0ToKpipipi.hh"
26#include "EvtGenBase/EvtComplex.hh"
27#include "EvtGenBase/EvtDecayTable.hh"
28#include <stdlib.h>
29
30EvtD0ToKpipipi::~EvtD0ToKpipipi() {}
31
32void EvtD0ToKpipipi::getName(std::string& model_name){
33 model_name="D0ToKpipipi";
34}
35
36EvtDecayBase* EvtD0ToKpipipi::clone(){
37 return new EvtD0ToKpipipi;
38}
39
40void EvtD0ToKpipipi::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 << "EvtD0ToKpipipi ==> Initialization !" << std::endl;
51
52 width[0] = 0.09;
53 width[1] = 0.044183653178315;
54 width[2] = 0.541879469380012;
55 width[3] = 0.148423336450619;
56 mass[0] = 1.272;
57 mass[1] = 0.894781734682169;
58 mass[2] = 1.3622013558915;
59 mass[3] = 0.779143408171384;
60
61 phi[0] = 2.34794687054858;
62 rho[0] = 0.0759345115620669;
63 phi[1] = -2.24641399153466;
64 rho[1] = 0.0383327604903577;
65 phi[2] = 2.48955684856045;
66 rho[2] = 0.0931445480476023;
67
68 phi[3] = 0;
69 rho[3] = 1;
70
71 phi[4] = -2.10558220063012;
72 rho[4] = 0.347041869435286;
73 phi[5] = 1.47445088061872;
74 phi[6] = 3.00243265559304;
75 rho[5] = 0.00965088341753795;
76 rho[6] = 0.120536507325731;
77 phi[7] = -2.45477499325158;
78 rho[7] = 0.101419048440676;
79 phi[8] = -1.35809992343491;
80 rho[8] = 4.28149643321317;
81 phi[9] = -2.45149221243198;
82 rho[9] = 0.339492272598394;
83 phi[10] = -0.17419389225461;
84 rho[10] = -0.143619437541254;
85
86 phi[11] = -2.08744386934208;
87 rho[11] = 0.296286583716349;
88 phi[12] = 0.;
89 rho[12] = 0.;
90 phi[13] = -0.432190571560873;
91 rho[13] = 0.657344690733276;
92 phi[14] = -1.39790294886865;
93 rho[14] = 1.71208007006123;
94 phi[15] = 1.58945300476228;
95 rho[15] = 3.58248347683687;
96 phi[16] = 2.58249107256307;
97 rho[16] = -1.10728829503506;
98 phi[17] = -0.163623135170955;
99 rho[17] = 1.70863070178363;
100 phi[18] = -0.134699023080211;
101 rho[18] = 0.567531283682344;
102 phi[19] = -2.12670610368279;
103 rho[19] = 0.276571752504914;
104 phi[20] = -1.3352622107357;
105 rho[20] = 0.416634203151278;
106
107 phi[21] = -2.91571684221842;
108 rho[21] = 0.423062298489176;
109 phi[22] = 2.4544220004327;
110 rho[22] = 1.4017194038459;
111 phi[23] = -2.23388390670423;
112 rho[23] = 4.11110400629068;
113
114 for (int i=0; i<24; i++) {
115 cout << i << "rho,phi = " << rho[i] << ", "<< phi[i] << endl;
116 }
117
118 mD = 1.86486;
119 rRes = 3.0;
120 rD = 5.0;
121 metap = 0.95778;
122 mkstr = 0.89594;
123 mk0 = 0.497614;
124 mass_Kaon = 0.49368;
125 mass_Pion = 0.13957;
126 mass_Pi0 = 0.1349766;
127 math_pi = 3.1415926;
128
129 pi = 3.1415926;
130 mpi = 0.13957;
131 g1 = 0.5468;
132 g2 = 0.23;
133
134 int GG[4][4] = { {1,0,0,0}, {0,-1,0,0}, {0,0,-1,0}, {0,0,0,-1} };
135 int EE[4][4][4][4] =
136 { { {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
137 {{0,0,0,0}, {0,0,0,0}, {0,0,0,1}, {0,0,-1,0}},
138 {{0,0,0,0}, {0,0,0,-1}, {0,0,0,0}, {0,1,0,0} },
139 {{0,0,0,0}, {0,0,1,0}, {0,-1,0,0}, {0,0,0,0} } },
140 { {{0,0,0,0}, {0,0,0,0}, {0,0,0,-1}, {0,0,1,0} },
141 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
142 {{0,0,0,1}, {0,0,0,0}, {0,0,0,0}, {-1,0,0,0}},
143 {{0,0,-1,0}, {0,0,0,0}, {1,0,0,0}, {0,0,0,0} } },
144 { {{0,0,0,0}, {0,0,0,1}, {0,0,0,0}, {0,-1,0,0}},
145 {{0,0,0,-1}, {0,0,0,0}, {0,0,0,0}, {1,0,0,0} },
146 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} },
147 {{0,1,0,0}, {-1,0,0,0}, {0,0,0,0}, {0,0,0,0} } },
148 { {{0,0,0,0}, {0,0,-1,0}, {0,1,0,0}, {0,0,0,0} },
149 {{0,0,1,0}, {0,0,0,0}, {-1,0,0,0}, {0,0,0,0} },
150 {{0,-1,0,0}, {1,0,0,0}, {0,0,0,0}, {0,0,0,0} },
151 {{0,0,0,0}, {0,0,0,0}, {0,0,0,0}, {0,0,0,0} } } };
152 for (int i=0; i<4; i++) {
153 for (int j=0; j<4; j++) {
154 G[i][j] = GG[i][j];
155 for (int k=0; k<4; k++) {
156 for (int l=0; l<4; l++) {
157 E[i][j][k][l] = EE[i][j][k][l];
158 }
159 }
160 }
161 }
162}
163
164void EvtD0ToKpipipi::initProbMax() {
165 setProbMax(720.0);
166}
167
168void EvtD0ToKpipipi::decay( EvtParticle *p){
169/*
170 double maxprob = 0.0;
171 for(int ir=0;ir<=60000000;ir++){
172 p->initializePhaseSpace(getNDaug(),getDaugs());
173 EvtVector4R Km0 = p->getDaug(0)->getP4();
174 EvtVector4R pi1 = p->getDaug(1)->getP4();
175 EvtVector4R pi2 = p->getDaug(2)->getP4();
176 EvtVector4R pi3 = p->getDaug(3)->getP4();
177 double Km[4],Pip1[4],Pip2[4],Pim[4];
178 Km[0] = Km0.get(0); Pip1[0] = pi1.get(0); Pip2[0] = pi2.get(0); Pim[0] = pi3.get(0);
179 Km[1] = Km0.get(1); Pip1[1] = pi1.get(1); Pip2[1] = pi2.get(1); Pim[1] = pi3.get(1);
180 Km[2] = Km0.get(2); Pip1[2] = pi1.get(2); Pip2[2] = pi2.get(2); Pim[2] = pi3.get(2);
181 Km[3] = Km0.get(3); Pip1[3] = pi1.get(3); Pip2[3] = pi2.get(3); Pim[3] = pi3.get(3);
182 double Prob = calPDF(Km, Pip1, Pip2, Pim);
183 if(Prob>maxprob) {
184 maxprob=Prob;
185 std::cout << "Max PDF = " << ir << " prob= " << Prob << std::endl;
186 }
187 }
188 std::cout << "Max!!!!!!!!!!! " << maxprob<< std::endl;
189*/
190 p->initializePhaseSpace(getNDaug(),getDaugs());
191 EvtVector4R Km0 = p->getDaug(0)->getP4();
192 EvtVector4R pi1 = p->getDaug(1)->getP4();
193 EvtVector4R pi2 = p->getDaug(2)->getP4();
194 EvtVector4R pi3 = p->getDaug(3)->getP4();
195
196 double Km[4],Pip1[4],Pip2[4],Pim[4];
197 Km[0] = Km0.get(0); Pip1[0] = pi1.get(0); Pip2[0] = pi2.get(0); Pim[0] = pi3.get(0);
198 Km[1] = Km0.get(1); Pip1[1] = pi1.get(1); Pip2[1] = pi2.get(1); Pim[1] = pi3.get(1);
199 Km[2] = Km0.get(2); Pip1[2] = pi1.get(2); Pip2[2] = pi2.get(2); Pim[2] = pi3.get(2);
200 Km[3] = Km0.get(3); Pip1[3] = pi1.get(3); Pip2[3] = pi2.get(3); Pim[3] = pi3.get(3);
201 double prob = calPDF(Km, Pip1, Pip2, Pim);
202 setProb(prob);
203 return;
204}
205
206double EvtD0ToKpipipi::calPDF(double Km[], double Pip1[], double Pip2[], double Pim[]) {
207 Km[0] = sqrt(mass_Kaon*mass_Kaon + Km[1]*Km[1] + Km[2]*Km[2] + Km[3]*Km[3]);
208 Pip1[0] = sqrt(mass_Pion*mass_Pion + Pip1[1]*Pip1[1] + Pip1[2]*Pip1[2] + Pip1[3]*Pip1[3]);
209 Pip2[0] = sqrt(mass_Pion*mass_Pion + Pip2[1]*Pip2[1] + Pip2[2]*Pip2[2] + Pip2[3]*Pip2[3]);
210 Pim[0] = sqrt(mass_Pion*mass_Pion + Pim[1]*Pim[1] + Pim[2]*Pim[2] + Pim[3]*Pim[3]);
211
212 EvtComplex PDF[24];
213 int g[3];
214 g[0] = 1; g[1] = 1;
215
216 //-----------D->K*rho--------
217 g[2] = 0;
218 PDF[0] = D2VV(Km,Pip1,Pip2,Pim,g) + D2VV(Km,Pip2,Pip1,Pim,g);
219 g[2] = 1;
220 PDF[1] = D2VV(Km,Pip1,Pip2,Pim,g) + D2VV(Km,Pip2,Pip1,Pim,g);
221 g[2] = 2;
222 PDF[2] = D2VV(Km,Pip1,Pip2,Pim,g) + D2VV(Km,Pip2,Pip1,Pim,g);
223 //-----------D->K*rho finish--
224 //----------D->a1K------------
225 g[2] = 0;
226 PDF[3] = D2AP_A2VP(Km,Pip2,Pip1,Pim,g,2) + D2AP_A2VP(Km,Pip1,Pip2,Pim,g,2);
227 g[2] = 2;
228 PDF[4] = D2AP_A2VP(Km,Pip2,Pip1,Pim,g,2) + D2AP_A2VP(Km,Pip1,Pip2,Pim,g,2);
229 //----------D->a1K finish-----
230 //--D->K1Pi--K1->K*Pi---------
231 g[2] = 0;
232 PDF[5] = D2AP_A2VP(Pip2,Pim,Km,Pip1,g,1) + D2AP_A2VP(Pip1,Pim,Km,Pip2,g,1);
233 g[2] = 2;
234 PDF[6] = D2AP_A2VP(Pip2,Pim,Km,Pip1,g,1) + D2AP_A2VP(Pip1,Pim,Km,Pip2,g,1);
235 //--D->K1Pi--K1->K*Pi-finish--
236 //--D->K1Pi--K1->rhoK---------
237 g[2] = 0;
238 PDF[7] = D2AP_A2VP(Pip2,Km,Pip1,Pim,g,3) + D2AP_A2VP(Pip1,Km,Pip2,Pim,g,3);
239 //--D->K1Pi--K1->rhoK-finish--
240 //--D->K1Pi--K1->K*0(1430)Pi--
241 PDF[8] = D2AP_A2SP(Pip2,Pim,Km,Pip1,1) + D2AP_A2SP(Pip1,Pim,Km,Pip2,1);
242 //--------res finish----------
243 //----------------non-res------------------
244 //--------KPiSrho------------------
245 PDF[9] = D2VS(Pip1,Pim,Km,Pip2,1,2) + D2VS(Pip2,Pim,Km,Pip1,1,2);
246 //--------K*PiPiS-----------------
247 PDF[10] = D2VS(Km,Pip1,Pip2,Pim,1,1) + D2VS(Km,Pip2,Pip1,Pim,1,1);
248 //--------K*PiP-------------------
249 PDF[11] = D2PP_P2VP(Pip2,Pim,Km,Pip1,1) + D2PP_P2VP(Pip1,Pim,Km,Pip2,1);
250 //--------rhoKA--------------------
251 g[0] = 1; g[1] = 0;
252 g[2] = 0;
253 PDF[12] = 0;
254 g[2] = 2;
255 PDF[13] = D2AP_A2VP(Pip2,Km,Pip1,Pim,g,3) + D2AP_A2VP(Pip1,Km,Pip2,Pim,g,3);
256 //-------PHSP-----------------------
257 PDF[14] = PHSP(Km,Pip1) + PHSP(Km,Pip2);
258 //------KPiVPiPiV-----------------------
259 g[0] = 0; g[1] = 0; g[2] = 0;
260 PDF[15] = D2VV(Km,Pip1,Pip2,Pim,g) + D2VV(Km,Pip2,Pip1,Pim,g);
261 //------KPiVPiPiS----------------------
262 PDF[16] = D2VS(Km,Pip1,Pip2,Pim,0,1) + D2VS(Km,Pip2,Pip1,Pim,0,1);
263 //------KPiSPiPiV----------------------
264 PDF[17] = D2VS(Pip1,Pim,Km,Pip2,0,2) + D2VS(Pip2,Pim,Km,Pip1,0,2);
265 //-------------------------------------
266 PDF[18] = D2PP_P2VP(Pip2,Km,Pip1,Pim,2) + D2PP_P2VP(Pip1,Km,Pip2,Pim,2);
267 //-------------------------------------
268 PDF[19] = D2VP_V2VP(Pip2,Pim,Km,Pip1,1) + D2VP_V2VP(Pip1,Pim,Km,Pip2,1);
269 //-------------------------------------
270 PDF[20] = D2VP_V2VP(Pip2,Km,Pip1,Pim,2) + D2VP_V2VP(Pip1,Km,Pip2,Pim,2);
271 //-------------------------------------
272 PDF[21] = D2TS(Km,Pip1,Pip2,Pim,1) + D2TS(Km,Pip2,Pip1,Pim,1);
273 PDF[22] = D2TS(Pip1,Pim,Km,Pip2,2) + D2TS(Pip2,Pim,Km,Pip1,2);
274 PDF[23] = D2AP_A2SP(Km,Pip2,Pip1,Pim,2) + D2AP_A2SP(Km,Pip1,Pip2,Pim,2);
275//------------------------------------------
276 EvtComplex cof;
277 EvtComplex pdf, module;
278 pdf = EvtComplex(0,0);
279 for(int i=0; i!=24; i++){
280 cof = EvtComplex(rho[i]*cos(phi[i]),rho[i]*sin(phi[i]));
281 //cout<<i << " " << (cof*PDF[i])<<" : "<<cof<<" "<<PDF[i]<<endl;
282 pdf = pdf + cof*PDF[i];
283 }
284 module = conj(pdf)*pdf;
285 double value;
286 value = real(module);
287 return (value <= 0) ? 1e-20 : value;
288}
289
290EvtComplex EvtD0ToKpipipi::KPiSFormfactor(double sa, double sb, double sc, double r)
291{
292 double m1430 = 1.463;
293 double sa0 = m1430*m1430;
294 double w1430 = 0.233;
295 double q0 = (sa0+sb-sc)*(sa0+sb-sc)/(4*sa0)-sb;
296 if(q0<0) q0 = 1e-16;
297 double qs = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
298 double q = sqrt(qs);
299 double width = w1430*q*m1430/sqrt(sa*q0);
300 double temp_R = atan(m1430*width/(sa0-sa));
301 if(temp_R<0) temp_R += math_pi;
302 double deltaR = -5.31 + temp_R;
303 double temp_F = atan(2*1.07*q/(2+(-1.8)*1.07*qs));
304 if(temp_F<0) temp_F += math_pi;
305 double deltaF = 2.33 + temp_F;
306 EvtComplex cR(cos(deltaR), sin(deltaR));
307 EvtComplex cF(cos(deltaF), sin(deltaF));
308 EvtComplex amp = 0.8*sin(deltaF)*cF + sin(deltaR)*cR*cF*cF;
309 return amp;
310}
311
312EvtComplex EvtD0ToKpipipi::D2VV(double P1[], double P2[], double P3[], double P4[], int g[])
313{
314 double t1V1[4], t1V2[4], t1D[4], t2D[4][4];
315 double temp_PDF = 0;
316 EvtComplex amp_PDF(0,0);
317 EvtComplex pro[2];
318 //---------use g[0],g[1] to define res or non-res, g[2] represent S, P or D wave
319 double sa[3], sb[3], sc[3], B[3];
320 double pV1[4], pV2[4], pD[4];
321 for(int i=0; i!=4; i++){
322 pV1[i] = P1[i] + P2[i];
323 pV2[i] = P3[i] + P4[i];
324 pD[i] = pV1[i] + pV2[i];
325 }
326 sa[0] = dot(pV1,pV1);
327 sb[0] = dot(P1,P1);
328 sc[0] = dot(P2,P2);
329 sa[1] = dot(pV2,pV2);
330 sb[1] = dot(P3,P3);
331 sc[1] = dot(P4,P4);
332 sa[2] = dot(pD,pD);
333 sb[2] = sa[0];
334 sc[2] = sa[1];
335 if(g[1] == 1){
336 pro[1] = propagatorGS(mass[3],width[3],sa[1],sb[1],sc[1],rRes,1);
337 }
338 if(g[0] == 1){
339 pro[0] = propagatorRBW(mass[1],width[1],sa[0],sb[0],sc[0],rRes,1);
340 }
341 if(g[0] == 0) pro[0] = 1;
342 if(g[1] == 0) pro[1] = 1;
343 B[0] = barrier(1,sa[0],sb[0],sc[0],rRes);
344 B[1] = barrier(1,sa[1],sb[1],sc[1],rRes);
345 calt1(P1,P2,t1V1);
346 calt1(P3,P4,t1V2);
347 if(g[2] == 0){
348 for(int i=0; i!=4; i++){
349 temp_PDF += (G[i][i])*t1V1[i]*t1V2[i];
350 }
351 B[2] = 1;
352 }
353 if(g[2] == 1){
354 calt1(pV1,pV2,t1D);
355 for(int i=0; i!=4; i++){
356 for(int j=0; j!=4; j++){
357 for(int k=0; k!=4; k++){
358 for(int l=0; l!=4; l++){
359 temp_PDF += E[i][j][k][l]*pD[i]*t1D[j]*t1V1[k]*t1V2[l]*(G[i][i])*(G[j][j])*(G[l][l])*(G[k][k]);
360 }
361 }
362 }
363 }
364 B[2] = barrier(1,sa[2],sb[2],sc[2],rD);
365 }
366 if(g[2] == 2){
367 calt2(pV1,pV2,t2D);
368 for(int i=0; i!=4; i++){
369 for(int j=0; j!=4; j++){
370 temp_PDF += t2D[i][j]*t1V1[i]*t1V2[j]*(G[i][i])*(G[j][j]);
371 }
372 }
373 B[2] = barrier(2,sa[2],sb[2],sc[2],rD);
374 }
375 amp_PDF = temp_PDF*B[0]*B[1]*B[2]*pro[0]*pro[1];
376 return amp_PDF;
377}
378
379EvtComplex EvtD0ToKpipipi::D2AP_A2VP(double P1[], double P2[], double P3[], double P4[], int g[], int flag)
380{
381 //flag = 1, V = K*, A = K1(1270); flag = 2, V = rho, A = a1(1260);
382 //flag = 3, V = rho, A = K1(1270);
383 double temp_PDF = 0;
384 EvtComplex amp_PDF(0,0);
385 EvtComplex pro[2];
386 double t1V[4], t1D[4], t2A[4][4];
387 double sa[3], sb[3], sc[3], B[3];
388 double pV[4],pA[4],pD[4];
389 for(int i=0; i!=4; i++){
390 pV[i] = P3[i] + P4[i];
391 pA[i] = pV[i] + P2[i];
392 pD[i] = pA[i] + P1[i];
393 }
394 sa[0] = dot(pV,pV);
395 sb[0] = dot(P3,P3);
396 sc[0] = dot(P4,P4);
397 sa[1] = dot(pA,pA);
398 sb[1] = sa[0];
399 sc[1] = dot(P2,P2);
400 sa[2] = dot(pD,pD);
401 sb[2] = sa[1];
402 sc[2] = dot(P1,P1);
403 if(g[0] == 1){
404 if(flag == 1) pro[0] = propagatorRBW(mass[1],width[1],sa[0],sb[0],sc[0],rRes,1);
405 if(flag == 2 || flag == 3) pro[0] = propagatorGS(mass[3],width[3],sa[0],sb[0],sc[0],rRes,1);
406 }
407 if(g[1] == 1){
408 if(flag == 1) pro[1] = propogator(mass[0],width[0],sa[1]);
409 if(flag == 2) pro[1] = propagatorRBW(mass[2],width[2],sa[1],sb[1],sc[1],rRes,g[2]);
410 if(flag == 3) pro[1] = propogator(mass[0],width[0],sa[1]);
411 }
412 if(g[0] == 0) pro[0] = 1;
413 if(g[1] == 0) pro[1] = 1;
414 B[0] = barrier(1,sa[0],sb[0],sc[0],rRes);
415 B[2] = barrier(1,sa[2],sb[2],sc[2],rD);
416 calt1(P3,P4,t1V);
417 calt1(pA,P1,t1D);
418 if(g[2] == 0){
419 for(int i=0; i!=4; i++){
420 for(int j=0; j!=4; j++){
421 temp_PDF += t1D[i]*(G[i][j] - pA[i]*pA[j]/sa[1])*t1V[j]*(G[i][i])*(G[j][j]);
422 }
423 }
424 B[1] = 1;
425 }
426 if(g[2] == 2){
427 calt2(pV,P2,t2A);
428 for(int i=0; i!=4; i++){
429 for(int j=0; j!=4; j++){
430 temp_PDF += t1D[i]*t2A[i][j]*t1V[j]*(G[i][i])*(G[j][j]);
431 }
432 }
433 B[1] = barrier(2,sa[1],sb[1],sc[1],rRes);
434 }
435 amp_PDF = temp_PDF*B[0]*B[1]*B[2]*pro[0]*pro[1];
436 return amp_PDF;
437}
438EvtComplex EvtD0ToKpipipi::D2AP_A2SP(double P1[], double P2[], double P3[], double P4[], int flag)
439{
440 //flag = 1, S = K*; flag = 2, S = rho
441 double temp_PDF = 0;
442 EvtComplex amp_PDF(0,0);
443 EvtComplex pro;
444 double sa[3], sb[3], sc[3], B[3];
445 double t1D[4], t1A[4];
446 double pS[4], pA[4], pD[4];
447 for(int i=0; i!=4; i++){
448 pS[i] = P3[i] + P4[i];
449 pA[i] = pS[i] + P2[i];
450 pD[i] = pA[i] + P1[i];
451 }
452 sa[0] = dot(pS,pS);
453 sb[0] = dot(P3,P3);
454 sc[0] = dot(P4,P4);
455 sa[1] = dot(pA,pA);
456 sb[1] = sa[0];
457 sc[1] = dot(P2,P2);
458 sa[2] = dot(pD,pD);
459 sb[2] = sa[1];
460 sc[2] = dot(P1,P1);
461 B[1] = barrier(1,sa[1],sb[1],sc[1],rRes);
462 B[2] = barrier(1,sa[2],sb[2],sc[2],rD);
463 calt1(pA,P1,t1D);
464 calt1(pS,P2,t1A);
465 for(int i=0; i!=4; i++){
466 temp_PDF += t1D[i]*t1A[i]*(G[i][i]);
467 }
468 amp_PDF = temp_PDF*B[1]*B[2];
469 if(flag == 1) amp_PDF *= KPiSFormfactor(sa[0],sb[0],sc[0],rRes);
470 return amp_PDF;
471}
472EvtComplex EvtD0ToKpipipi::D2PP_P2VP(double P1[], double P2[], double P3[], double P4[], int flag)
473{
474 //flag = 1, V = K*; flag = 2, V = rho
475 double temp_PDF = 0;
476 EvtComplex amp_PDF(0,0);
477 EvtComplex pro;
478 double sa[3], sb[3], sc[3], B[3];
479 double t1V[4];
480 double pV[4], pP[4], pD[4];
481 for(int i=0; i!=4; i++){
482 pV[i] = P3[i] + P4[i];
483 pP[i] = pV[i] + P2[i];
484 pD[i] = pP[i] + P1[i];
485 }
486 sa[0] = dot(pV,pV);
487 sb[0] = dot(P3,P3);
488 sc[0] = dot(P4,P4);
489 sa[1] = dot(pP,pP);
490 sb[1] = sa[0];
491 sc[1] = dot(P2,P2);
492 sa[2] = dot(pD,pD);
493 sb[2] = sa[1];
494 sc[2] = dot(P1,P1);
495 B[0] = barrier(1,sa[0],sb[0],sc[0],rRes);
496 B[1] = barrier(1,sa[1],sb[1],sc[1],rRes);
497 if(flag == 1) pro = propagatorRBW(mass[1],width[1],sa[0],sb[0],sc[0],rRes,1);
498 if(flag == 2) pro = propagatorGS(mass[3],width[3],sa[0],sb[0],sc[0],rRes,1);
499 calt1(P3,P4,t1V);
500 for(int i=0; i!=4; i++){
501 temp_PDF += P2[i]*t1V[i]*(G[i][i]);
502 }
503 amp_PDF = temp_PDF*B[0]*B[1]*pro;
504 return amp_PDF;
505}
506EvtComplex EvtD0ToKpipipi::D2VP_V2VP(double P1[], double P2[], double P3[], double P4[], int flag)
507{
508 //flag = 1, (K*Pi)V; flag = 2, (rhoK)V
509 double temp_PDF = 0;
510 EvtComplex amp_PDF(0,0);
511 EvtComplex pro;
512 double sa[3], sb[3], sc[3], B[3];
513 double pV1[4], pV2[4], qV1[4], qV2[4], pD[4];
514 for(int i=0; i!=4; i++){
515 pV2[i] = P3[i] + P4[i];
516 qV2[i] = P3[i] - P4[i];
517 pV1[i] = pV2[i] + P2[i];
518 qV1[i] = pV2[i] - P2[i];
519 pD[i] = pV1[i] + P1[i];
520 }
521 for(int i=0; i!=4; i++){
522 for(int j=0; j!=4; j++){
523 for(int k=0; k!=4; k++){
524 for(int l=0; l!=4; l++){
525 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]);
526 }
527 }
528 }
529 }
530 sa[0] = dot(pV2,pV2);
531 sb[0] = dot(P3,P3);
532 sc[0] = dot(P4,P4);
533 sa[1] = dot(pV1,pV1);
534 sb[1] = sa[0];
535 sc[1] = dot(P2,P2);
536 sa[2] = dot(pD,pD);
537 sb[2] = sa[1];
538 sc[2] = dot(P1,P1);
539 if(flag == 1) pro = propagatorRBW(mass[1],width[1],sa[0],sb[0],sc[0],rRes,1);
540 if(flag == 2) pro = propagatorGS(mass[3],width[3],sa[0],sb[0],sc[0],rRes,1);
541 B[0] = barrier(1,sa[0],sb[0],sc[0],rRes);
542 B[1] = barrier(1,sa[1],sb[1],sc[1],rRes);
543 B[2] = barrier(1,sa[2],sb[2],sc[2],rD);
544 amp_PDF = temp_PDF*B[0]*B[1]*B[2]*pro;
545 return amp_PDF;
546}
547EvtComplex EvtD0ToKpipipi::D2VS(double P1[], double P2[], double P3[], double P4[], int g, int flag)
548{
549 //flag = 1, V = K*; flag = 2, V = rho
550 double temp_PDF = 0;
551 EvtComplex amp_PDF(0,0);
552 EvtComplex pro;
553 double sa[3], sb[3], sc[3], B[3];
554 double t1D[4], t1V[4];
555 double pS[4], pV[4], pD[4];
556 for(int i=0; i!=4; i++){
557 pS[i] = P3[i] + P4[i];
558 pV[i] = P1[i] + P2[i];
559 pD[i] = pS[i] + pV[i];
560 }
561 sa[0] = dot(pS,pS);
562 sb[0] = dot(P3,P3);
563 sc[0] = dot(P4,P4);
564 sa[1] = dot(pV,pV);
565 sb[1] = dot(P1,P1);
566 sc[1] = dot(P2,P2);
567 sa[2] = dot(pD,pD);
568 sb[2] = sa[0];
569 sc[2] = sa[1];
570 if(g == 1){
571 if(flag == 2) pro = propagatorGS(mass[3],width[3],sa[1],sb[1],sc[1],rRes,1);
572 if(flag == 1) pro = propagatorRBW(mass[1],width[1],sa[1],sb[1],sc[1],rRes,1);
573 }
574 if(g == 0) pro = 1;
575 B[1] = barrier(1,sa[1],sb[1],sc[1],rRes);
576 B[2] = barrier(1,sa[2],sb[2],sc[2],rD);
577 calt1(P1,P2,t1V);
578 calt1(pS,pV,t1D);
579 for(int i=0; i!=4; i++){
580 temp_PDF += G[i][i]*t1D[i]*t1V[i];
581 }
582 amp_PDF = temp_PDF*B[1]*B[2]*pro;
583 if(flag == 2) amp_PDF *= KPiSFormfactor(sa[0],sb[0],sc[0],rRes);
584 return amp_PDF;
585}
586EvtComplex EvtD0ToKpipipi::D2TS(double P1[], double P2[], double P3[], double P4[], int flag)
587{
588 //flag = 1, T = K*; flag = 2, T = rho
589 double temp_PDF = 0;
590 EvtComplex amp_PDF(0,0);
591 double sa[3], sb[3], sc[3], B[3];
592 double t2D[4][4], t2T[4][4];
593 double pS[4], pT[4], pD[4];
594 for(int i=0; i!=4; i++){
595 pS[i] = P3[i] + P4[i];
596 pT[i] = P1[i] + P2[i];
597 pD[i] = pT[i] + pS[i];
598 }
599 sa[0] = dot(pT,pT);
600 sb[0] = dot(P1,P1);
601 sc[0] = dot(P2,P2);
602 sa[1] = dot(pS,pS);
603 sb[1] = dot(P3,P3);
604 sc[1] = dot(P4,P4);
605 sa[2] = dot(pD,pD);
606 sb[2] = sa[0];
607 sc[2] = sa[1];
608 B[0] = barrier(2,sa[0],sb[0],sc[0],rRes);
609 B[2] = barrier(2,sa[2],sb[2],sc[2],rD);
610 calt2(P1,P2,t2T);
611 calt2(pT,pS,t2D);
612 for(int i=0; i!=4; i++){
613 for(int j=0; j!=4; j++){
614 temp_PDF += t2D[i][j]*t2T[j][i]*(G[i][i])*(G[j][j]);
615 }
616 }
617 amp_PDF = temp_PDF*B[0]*B[2];
618 if(flag == 2) amp_PDF *= KPiSFormfactor(sa[1],sb[1],sc[1],rRes);
619 return amp_PDF;
620}
621EvtComplex EvtD0ToKpipipi::PHSP(double Km[], double Pip[])
622{
623 EvtComplex amp_PDF(0,0);
624 double sa,sb,sc;
625 double KPi[4];
626 for(int i=0; i!=4; i++){
627 KPi[i] = Km[i] + Pip[i];
628 }
629 sa = dot(KPi,KPi);
630 sb = dot(Km,Km);
631 sc = dot(Pip,Pip);
632 amp_PDF = KPiSFormfactor(sa,sb,sc,rRes);
633 return amp_PDF;
634}
635double EvtD0ToKpipipi::calDalEva(double P1[], double P2[], double P3[])
636{
637 //Ks Pi0 Pi0 exchange Pi0 in this case
638 //CLEOc model
639 EvtComplex PDF[7], cof, pdf, module;
640 //Ks Pi0 Pi0
641 double P[3][4];
642 for(int i=0; i<4; i++){
643 P[0][i] = P1[i];
644 P[1][i] = P2[i];
645 P[2][i] = P3[i];
646 }
647 PDF[0] = Spin_factor(P[1],P[2],P[0],0) + Spin_factor(P[2],P[1],P[0],0);
648 PDF[1] = Spin_factor(P[1],P[2],P[0],10) + Spin_factor(P[2],P[1],P[0],10);
649 PDF[2] = Spin_factor(P[1],P[2],P[0],100) + Spin_factor(P[2],P[1],P[0],100);
650 //-----------------
651 PDF[3] = Spin_factor(P[0],P[1],P[2],1) + Spin_factor(P[0],P[2],P[1],1);
652 PDF[4] = Spin_factor(P[0],P[1],P[2],11) + Spin_factor(P[0],P[2],P[1],11);
653 //-----------------
654 PDF[5] = Spin_factor(P[1],P[2],P[0],2) + Spin_factor(P[2],P[1],P[0],2);
655 PDF[6] = Spin_factor(P[0],P[1],P[2],12) + Spin_factor(P[0],P[2],P[1],12);
656
657 rho[0] = 4.2551e-01; phi[0] = 3.7370e+00;
658 rho[1] = -1.0086e+00; phi[1] = 6.9153e+00;
659 rho[2] = 6.5290e-01; phi[2] = 2.8041e+00;
660 rho[3] = 1.0; phi[3] = 0.0;
661 rho[4] = -2.5079e-01; phi[4] = 6.3740e-01;
662 rho[5] = -4.1640e-01; phi[5] = 7.7850e+00;
663 rho[6] = 2.1281e-01; phi[6] = 5.3474e+00;
664
665 pdf = EvtComplex(0,0);
666 for(int i=0; i<7; i++){
667 cof = EvtComplex(rho[i]*cos(phi[i]),rho[i]*sin(phi[i]));
668 pdf += cof*PDF[i];
669 }
670 module = conj(pdf)*pdf;
671 double value;
672 value = real(module);
673 return (value <= 0) ? 1e-20 : value;
674}
675EvtComplex EvtD0ToKpipipi::Spin_factor(double P1[], double P2[], double P3[], int spin)
676{
677 // D-> R P3, R->P1 P2
678 double R[4], s[3], sp2, sD, B[2];
679 for(int i=0; i<4; i++){
680 R[i] = P1[i] + P2[i];
681 }
682 s[0] = dot(R,R);
683 s[1] = dot(P1, P1);
684 s[2] = dot(P2, P2);
685 sp2 = dot(P3,P3);
686 sD = mD*mD;
687 EvtComplex amp(0,0), prop;
688 prop = getProp(s, spin);
689 double tmp(0.);
690 if(spin%10 == 0){
691 amp = prop;
692 }
693 else if(spin%10 == 1){
694 tmp = 0;
695 double T1[4], t1[4];
696 calt1(R,P3,T1);
697 calt1(P1,P2,t1);
698 B[0] = barrier(1,s[0],s[1],s[2],3.0);
699 B[1] = barrier(1,sD,s[0],sp2,3.0);
700 for(int i=0; i<4; i++){
701 tmp += T1[i]*t1[i]*G[i][i];
702 }
703 amp = tmp*prop*B[0]*B[1];
704 }
705 else if(spin%10 ==2){
706 tmp = 0;
707 double T2[4][4], t2[4][4];
708 calt2(R,P3,T2);
709 calt2(P1,P2,t2);
710 B[0] = barrier(2,s[0],s[1],s[2],3.0);
711 B[1] = barrier(2,sD,s[0],sp2,3.0);
712 for(int i=0; i<4; i++){
713 for(int j=0; j<4; j++){
714 tmp += T2[i][j]*t2[j][i]*G[j][j]*G[i][i];
715 }
716 }
717 amp = tmp*prop*B[0]*B[1];
718 }
719 else{
720 cout<<"Only S, P, D wave allowed"<<endl;
721 }
722 return amp;
723}
724EvtComplex EvtD0ToKpipipi::getProp(double s[], int flag)
725{
726 //Ks Pi0 Pi0
727 //0 980
728 //10 f0_1370
729 //20 f0_1500
730 //2 f2_1270
731 //1 K*892
732 //12 K*2_1430
733 //11 K*1680
734 //100 f0(500)
735 EvtComplex prop(1,0);
736 double snpi, scpi, snk, sck;
737 snpi = mass_Pi0*mass_Pi0;
738 scpi = mass_Pion*mass_Pion;
739 sck = mass_Kaon*mass_Kaon;
740 snk = mk0*mk0;
741 if(flag == 0){
742 double mass, gpipi, gkk;
743 mass = 0.965;
744 gpipi = 0.406;
745 gkk = 2*gpipi;
746 EvtComplex ci(0,1);
747 EvtComplex rhokk, rhopipi;
748 rhokk = 0.5*(rhofactor(s[0],sck)+rhofactor(s[0],snk));
749 rhopipi = 1.0/3.0 * (2.0 * rhofactor(s[0],scpi) + rhofactor(s[0],snpi));
750 prop = 1.0/(mass*mass - s[0] - ci*(gpipi*gpipi*rhopipi+gkk*gkk*rhokk));
751 }
752 if(flag == 10){
753 prop = propagatorRBW(1.35,0.265,s[0],s[1],s[2],3.0,0);
754 }
755 if(flag == 20){
756 prop = propagatorRBW(1.505,0.109,s[0],s[1],s[2],3.0,0);
757 }
758 if(flag == 1){
759 prop = propagatorRBW(0.896,0.0503,s[0],s[1],s[2],3.0,1);
760 }
761 if(flag == 11){
762 prop = propagatorRBW(1.717,0.322,s[0],s[1],s[2],3.0,1);
763 }
764 if(flag == 2){
765 prop = propagatorRBW(1.2751,0.185,s[0],s[1],s[2],3.0,2);
766 }
767 if(flag == 12){
768 prop = propagatorRBW(1.4324,0.109,s[0],s[1],s[2],3.0,2);
769 }
770 if(flag == 100){
771 EvtComplex pole(0.470,-0.220);
772 prop = 1.0/(s[0] - pole*pole);
773 }
774 return prop;
775}
776EvtComplex EvtD0ToKpipipi::rhofactor(double sx, double sdau)
777{
778 EvtComplex one(1,0);
779 EvtComplex ci(0,1);
780 EvtComplex res;
781 double tmp;
782 tmp = 1 - 4*sdau/sx;
783 if(tmp>0) res = one*sqrt(tmp);
784 if(tmp<0) res = ci*sqrt(fabs(tmp));
785 return res;
786}
787EvtComplex EvtD0ToKpipipi::propogator(double mass, double width, double sx)const
788{
789 EvtComplex ci(0,1);
790 EvtComplex prop=1.0/(mass*mass-sx-ci*mass*width);
791 return prop;
792}
793double EvtD0ToKpipipi::wid(double mass, double sa, double sb, double sc, double r, int l)const
794{
795 double widm(0.), q(0.), q0(0.);
796 double sa0 = mass*mass;
797 double m = sqrt(sa);
798 q = Qabcs(sa,sb,sc);
799 q0 = Qabcs(sa0,sb,sc);
800 double z,z0;
801 z = q*r*r;
802 z0 = q0*r*r;
803 double t = q/q0;
804 double F(0.);
805 if(l == 0) F = 1;
806 if(l == 1) F = sqrt((1+z0)/(1+z));
807 if(l == 2) F = sqrt((9+3*z0+z0*z0)/(9+3*z+z*z));
808 widm = pow(t,l+0.5)*mass/m*F*F;
809 return widm;
810}
811double EvtD0ToKpipipi::h(double m, double q)const
812{
813 double h(0.);
814 h = 2/pi*q/m*log((m+2*q)/(2*mpi));
815 return h;
816}
817double EvtD0ToKpipipi::dh(double mass, double q0)const
818{
819 double dh = h(mass,q0)*(1.0/(8*q0*q0)-1.0/(2*mass*mass))+1.0/(2*pi*mass*mass);
820 return dh;
821}
822double EvtD0ToKpipipi::f(double mass, double sx, double q0, double q)const
823{
824 double m = sqrt(sx);
825 double f = mass*mass/(pow(q0,3))*(q*q*(h(m,q)-h(mass,q0))+(mass*mass-sx)*q0*q0*dh(mass,q0));
826 return f;
827}
828double EvtD0ToKpipipi::d(double mass, double q0)const
829{
830 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));
831 return d;
832}
833EvtComplex EvtD0ToKpipipi::propagatorRBW(double mass, double width, double sa, double sb, double sc, double r, int l)const
834{
835 EvtComplex ci(0,1);
836 EvtComplex prop=1.0/(mass*mass-sa-ci*mass*width*wid(mass,sa,sb,sc,r,l));
837 return prop;
838}
839EvtComplex EvtD0ToKpipipi::propagatorGS(double mass, double width, double sa, double sb, double sc, double r, int l)const
840{
841 EvtComplex ci(0,1);
842 double q = Qabcs(sa,sb,sc);
843 double sa0 = mass*mass;
844 double q0 = Qabcs(sa0,sb,sc);
845 q = sqrt(q);
846 q0 = sqrt(q0);
847 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));
848 return prop;
849}
850double EvtD0ToKpipipi::Flatte_rhoab(double sa, double sb, double sc)const
851{
852 double q = Qabcs(sa,sb,sc);
853 double rho = sqrt(q/sa);
854 return rho;
855}
856EvtComplex EvtD0ToKpipipi::propagatorFlatte(double mass, double width, double sx, double *sb, double *sc)const
857{
858 EvtComplex ci(0,1);
859 double rho1 = Flatte_rhoab(sx,sb[0],sc[0]);
860 double rho2 = Flatte_rhoab(sx,sb[1],sc[1]);
861 EvtComplex prop = 1.0/(mass*mass-sx-ci*(g1*g1*rho1+g2*g2*rho2));
862 return prop;
863}
864double EvtD0ToKpipipi::dot(double *a1, double *a2)const
865{
866 double dot = 0;
867 for(int i=0; i!=4; i++){
868 dot += a1[i]*a2[i]*G[i][i];
869 }
870 return dot;
871}
872double EvtD0ToKpipipi::Qabcs(double sa, double sb, double sc)const
873{
874 double Qabcs = (sa+sb-sc)*(sa+sb-sc)/(4*sa)-sb;
875 if(Qabcs < 0) Qabcs = 1e-16;
876 return Qabcs;
877}
878double EvtD0ToKpipipi::barrier(double l, double sa, double sb, double sc, double r)const
879{
880 double q = Qabcs(sa,sb,sc);
881 q = sqrt(q);
882 double z = q*r;
883 z = z*z;
884 double F = 1;
885 if(l > 2) F = 0;
886 if(l == 0)
887 F = 1;
888 if(l == 1){
889 F = sqrt((2*z)/(1+z));
890 }
891 if(l == 2){
892 F = sqrt((13*z*z)/(9+3*z+z*z));
893 }
894 return F;
895}
896void EvtD0ToKpipipi::calt1(double daug1[], double daug2[], double t1[]) const
897{
898 double p, pq;
899 double pa[4], qa[4];
900 for(int i=0; i!=4; i++){
901 pa[i] = daug1[i] + daug2[i];
902 qa[i] = daug1[i] - daug2[i];
903 }
904 p = dot(pa,pa);
905 pq = dot(pa,qa);
906 for(int i=0; i!=4; i++){
907 t1[i] = qa[i] - pq/p*pa[i];
908 }
909}
910void EvtD0ToKpipipi::calt2(double daug1[], double daug2[], double t2[][4]) const
911{
912 double p,r;
913 double pa[4], t1[4];
914 calt1(daug1,daug2,t1);
915 r = dot(t1,t1);
916 for(int i=0; i!=4; i++){
917 pa[i] = daug1[i] + daug2[i];
918 }
919 p = dot(pa,pa);
920 for(int i=0; i!=4; i++){
921 for(int j=0; j!=4; j++){
922 t2[i][j] = t1[i]*t1[j] - 1.0/3*r*(G[i][j]-pa[i]*pa[j]/p);
923 }
924 }
925}
sin
double sin(const BesAngle a)
Definition: BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: BesAngle.h:213
P
double P(RecMdcKalTrack *trk)
Definition: CalibEventSelect.cxx:87
mass
double mass
Definition: CosmicGenerator.cxx:138
f
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
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
conj
Evt3Rank3C conj(const Evt3Rank3C &t2)
Definition: bak-BesEvtGen-00-04-08/src/EvtGen/EvtGenBase/Evt3Rank3C.cc:175
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::R
complex_t R(double Q2, double M2, double G, double Mp2, double Mm2)
Definition: TUtil.h:27