BOSS 7.0.9
BESIII Offline Software System
Loading...
Searching...
No Matches
BesEvtGen-00-04-08/src/EvtGen/EvtGenModels/EvtDToKpienu.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+ -> K- pi+ e+ nu,
14// PHYSICAL REVIEW D 94, 032001 (2016)
15//
16// Modification history:
17// Liaoyuan Dong Jun.24, 2019 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/EvtDToKpienu.hh"
26#include "EvtGenBase/EvtComplex.hh"
27#include "EvtGenBase/EvtDecayTable.hh"
28#include <stdlib.h>
29
30EvtDToKpienu::~EvtDToKpienu() {}
31
32void EvtDToKpienu::getName(std::string& model_name){
33 model_name="DToKpienu";
34}
35
36EvtDecayBase* EvtDToKpienu::clone(){
37 return new EvtDToKpienu;
38}
39
40void EvtDToKpienu::init(){
41 checkNArg(0);
42 checkNDaug(4);
43 checkSpinParent(EvtSpinType::SCALAR);
44 checkSpinDaughter(0,EvtSpinType::SCALAR);
45 checkSpinDaughter(1,EvtSpinType::SCALAR);
46
47 std::cout << "EvtDToKpienu ==> Initialization !" << std::endl;
48 nAmps = 2;
49
50 rS = -11.57; // S-wave
51 rS1 = 0.08;
52 a_delta = 1.94;
53 b_delta = -0.81;
54 m0_1430_S = 1.425;
55 width0_1430_S = 0.270;
56 type[0] = 0;
57
58 mV = 1.81;
59 mA = 2.61;
60 V_0 = 1.411;
61 A1_0 = 1;
62 A2_0 = 0.788;
63
64 m0 = 0.8946; // P-wave K*
65 width0 = 0.04642;
66 rBW = 3.07;
67 rho = 1;
68 phi = 0;
69 type[1] = 1;
70
71 m0_1410 = 1.414; // P-wave K*(1410)
72 width0_1410 = 0.232;
73 rho_1410 = 0.1;
74 phi_1410 = 0.;
75 type[2] = 2;
76
77 TV_0 = 1; // D-wave K*2(1430)
78 T1_0 = 1;
79 T2_0 = 1;
80 m0_1430 = 1.4324;
81 width0_1430 = 0.109;
82 rho_1430 = 15;
83 phi_1430 = 0;
84 type[3] = 3;
85
86 mD = 1.86962;
87 mPi = 0.13957;
88 mK = 0.49368;
89 Pi = atan2(0.0,-1.0);
90 root2 = sqrt(2.);
91 root2d3 = sqrt(2./3);
92 root1d2 = sqrt(0.5);
93 root3d2 = sqrt(1.5);
94}
95
96void EvtDToKpienu::initProbMax() {
97 setProbMax(22750.0);
98}
99
100void EvtDToKpienu::decay( EvtParticle *p){
101/*
102 double maxprob = 0.0;
103 for(int ir=0;ir<=60000000;ir++){
104 p->initializePhaseSpace(getNDaug(),getDaugs());
105 EvtVector4R _K = p->getDaug(0)->getP4();
106 EvtVector4R _pi = p->getDaug(1)->getP4();
107 EvtVector4R _e = p->getDaug(2)->getP4();
108 EvtVector4R _nu = p->getDaug(3)->getP4();
109 int pid = EvtPDL::getStdHep(p->getDaug(0)->getId());
110 int charm;
111 if(pid == -321) charm = 1;
112 else charm = -1;
113 double m2, q2, cosV, cosL, chi;
114 KinVGen(_K, _pi, _e, _nu, charm, m2, q2, cosV, cosL, chi);
115 double _prob = calPDF(m2, q2, cosV, cosL, chi);
116 if(_prob>maxprob) {
117 maxprob=_prob;
118 std::cout << "Max PDF = " << ir << " charm= " << charm << " prob= " << _prob << std::endl;
119 }
120 }
121 std::cout << "Max!!!!!!!!!!! " << maxprob<< std::endl;
122*/
123 p->initializePhaseSpace(getNDaug(),getDaugs());
124 EvtVector4R K = p->getDaug(0)->getP4();
125 EvtVector4R pi = p->getDaug(1)->getP4();
126 EvtVector4R e = p->getDaug(2)->getP4();
127 EvtVector4R nu = p->getDaug(3)->getP4();
128
129 int pid = EvtPDL::getStdHep(p->getDaug(0)->getId());
130 int charm;
131 if(pid == -321) charm = 1;
132 else charm = -1;
133 double m2, q2, cosV, cosL, chi;
134 KinVGen(K, pi, e, nu, charm, m2, q2, cosV, cosL, chi);
135 double prob = calPDF(m2, q2, cosV, cosL, chi);
136 setProb(prob);
137 return;
138}
139
140void EvtDToKpienu::KinVGen(EvtVector4R vp4_K, EvtVector4R vp4_Pi, EvtVector4R vp4_Lep, EvtVector4R vp4_Nu, int charm, double& m2, double& q2, double& cosV, double& cosL, double& chi)
141{
142 EvtVector4R vp4_KPi = vp4_K + vp4_Pi;
143 EvtVector4R vp4_W = vp4_Lep + vp4_Nu;
144
145 m2 = vp4_KPi.mass2();
146 q2 = vp4_W.mass2();
147
148 EvtVector4R boost;
149 boost.set(vp4_KPi.get(0), -vp4_KPi.get(1), -vp4_KPi.get(2), -vp4_KPi.get(3));
150 EvtVector4R vp4_Kp = boostTo(vp4_K, boost);
151 cosV = vp4_Kp.dot(vp4_KPi)/(vp4_Kp.d3mag()*vp4_KPi.d3mag());
152
153 boost.set(vp4_W.get(0), -vp4_W.get(1), -vp4_W.get(2), -vp4_W.get(3));
154 EvtVector4R vp4_Lepp = boostTo(vp4_Lep, boost);
155 cosL = vp4_Lepp.dot(vp4_W)/(vp4_Lepp.d3mag()*vp4_W.d3mag());
156
157 EvtVector4R V = vp4_KPi/vp4_KPi.d3mag();
158 EvtVector4R C = vp4_Kp.cross(V);
159 C/=C.d3mag();
160 EvtVector4R D = vp4_Lepp.cross(V);
161 D/=D.d3mag();
162 double sinx = C.cross(V).dot(D);
163 double cosx = C.dot(D);
164 chi = sinx>0? acos(cosx): -acos(cosx);
165 if(charm==-1) chi=-chi;
166}
167
168double EvtDToKpienu::calPDF(double m2, double q2, double cosV, double cosL, double chi) {
169 double delta[5] = {0};
170 double amplitude[5] = {0};
171 double m = sqrt(m2);
172 double q = sqrt(q2);
173
174 //begin to calculate form factor
175 EvtComplex F10(0.0,0.0);
176 EvtComplex F11(0.0,0.0);
177 EvtComplex F21(0.0,0.0);
178 EvtComplex F31(0.0,0.0);
179 EvtComplex F12(0.0,0.0);
180 EvtComplex F22(0.0,0.0);
181 EvtComplex F32(0.0,0.0);
182
183 EvtComplex f10(0.0,0.0);
184 EvtComplex f11(0.0,0.0);
185 EvtComplex f21(0.0,0.0);
186 EvtComplex f31(0.0,0.0);
187 EvtComplex f12(0.0,0.0);
188 EvtComplex f22(0.0,0.0);
189 EvtComplex f32(0.0,0.0);
190 EvtComplex coef(0.0,0.0);
191 double amplitude_temp, delta_temp;
192
193 for(int index=0; index<nAmps; index++)
194 {
195 switch(type[index])
196 {
197 case 0: //calculate form factor of S wave
198 {
199 NRS(m, q, rS, rS1, a_delta, b_delta, mA, m0_1430_S, width0_1430_S, amplitude_temp, delta_temp, f10);
200 amplitude[index] = amplitude_temp;
201 delta[index] = delta_temp;
202 F10 = F10 + f10;
203 break;
204 }
205 case 1: //calculate form factor of P wave (K*)
206 {
207 ResonanceP(m, q, mV, mA, V_0, A1_0, A2_0, m0, width0, rBW, amplitude_temp, delta_temp, f11, f21, f31);
208 amplitude[index] = amplitude_temp;
209 delta[index] = delta_temp;
210 coef = getCoef(rho, phi);
211 F11 = F11+ coef*f11;
212 F21 = F21+ coef*f21;
213 F31 = F31+ coef*f31;
214 break;
215 }
216 case 2: //calculate form factor of P wave (K*(1410))
217 {
218
219 ResonanceP(m, q, mV, mA, V_0, A1_0, A2_0, m0_1410, width0_1410, rBW, amplitude_temp, delta_temp, f11, f21, f31);
220 amplitude[index] = amplitude_temp;
221 delta[index] = delta_temp;
222 coef = getCoef(rho_1410, phi_1410);
223 F11 = F11+ coef*f11;
224 F21 = F21+ coef*f21;
225 F31 = F31+ coef*f31;
226 break;
227 }
228 case 3: //calculate form factor of D wave
229 {
230 ResonanceD(m, q, mV, mA, TV_0, T1_0, T2_0, m0_1430, width0_1430, rBW, amplitude_temp, delta_temp, f12, f22, f32);
231 amplitude[index] = amplitude_temp;
232 delta[index] = delta_temp;
233 coef = getCoef(rho_1430, phi_1430);
234 F12 = F12+ coef*f12;
235 F22 = F22+ coef*f22;
236 F32 = F32+ coef*f32;
237 break;
238 }
239 default:
240 {
241 std::cout<<"No such form factor type!!!"<<std::endl;
242 break;
243 }
244 }
245 }
246
247 //begin to calculate pdf value
248 double I,I1,I2,I3,I4,I5,I6,I7,I8,I9;
249
250 double cosV2 = cosV*cosV;
251 double sinV = sqrt(1.0-cosV2);
252 double sinV2 = sinV*sinV;
253
254 EvtComplex F1 = F10 + F11*cosV + F12*(1.5*cosV2-0.5);
255 EvtComplex F2 = F21*root1d2 + F22*cosV*root3d2;
256 EvtComplex F3 = F31*root1d2 + F32*cosV*root3d2;
257
258 I1 = 0.25*( abs2(F1) + 1.5*sinV2*(abs2(F2) + abs2(F3) ) );
259 I2 = -0.25*( abs2(F1) - 0.5*sinV2*(abs2(F2) + abs2(F3) ) );
260 I3 = -0.25*( abs2(F2) - abs2(F3) )*sinV2;
261 I4 = real( conj(F1)*F2 )*sinV*0.5;
262 I5 = real( conj(F1)*F3 )*sinV;
263 I6 = real( conj(F2)*F3 )*sinV2;
264 I7 = imag( conj(F2)*F1 )*sinV;
265 I8 = imag( conj(F3)*F1 )*sinV*0.5;
266 I9 = imag( conj(F3)*F2 )*sinV2*(-0.5);
267
268 double coschi = cos(chi);
269 double sinchi = sin(chi);
270 double sin2chi = 2.0*sinchi*coschi;
271 double cos2chi = 1.0 - 2.0*sinchi*sinchi;
272
273 double sinL = sqrt(1.-cosL*cosL);
274 double sinL2 = sinL*sinL;
275 double sin2L = 2.0*sinL*cosL;
276 double cos2L = 1.0 - 2.0*sinL2;
277
278 I = I1 + I2*cos2L + I3*sinL2*cos2chi + I4*sin2L*coschi + I5*sinL*coschi + I6*cosL + I7*sinL*sinchi + I8*sin2L*sinchi + I9*sinL2*sin2chi;
279 return I;
280}
281
282void EvtDToKpienu::ResonanceP(double m, double q, double mV, double mA, double V_0, double A1_0, double A2_0, double m0, double width0, double rBW, double& amplitude, double& delta, EvtComplex& F11, EvtComplex& F21, EvtComplex& F31)
283{
284 double pKPi = getPStar(mD, m, q);
285 double mD2 = mD*mD;
286 double m2 = m*m;
287 double m02 = m0*m0;
288 double q2 = q*q;
289 double mV2 = mV*mV;
290 double mA2 = mA*mA;
291 double summDm = mD+m;
292 double V = V_0 /(1.0-q2/(mV2));
293 double A1 = A1_0/(1.0-q2/(mA2));
294 double A2 = A2_0/(1.0-q2/(mA2));
295 double A = summDm*A1;
296 double B = 2.0*mD*pKPi/summDm*V;
297
298 //construct the helicity form factor
299 double H0 = 0.5/(m*q)*((mD2-m2-q2)*summDm*A1-4.0*(mD2*pKPi*pKPi)/summDm*A2);
300 double Hp = A-B;
301 double Hm = A+B;
302
303 //calculate alpha
304 double B_Kstar = 2./3.; //B_Kstar = Br(Kstar(892)->k pi)
305 double pStar0 = getPStar(m0, mPi, mK);
306 double alpha = sqrt(3.*Pi*B_Kstar/(pStar0*width0));
307
308 //construct amplitudes of (non)resonance
309 double F = getF1(m, m0, mPi, mK, rBW);
310 double width = getWidth1(m, m0, mPi, mK, width0, rBW);
311
312 EvtComplex C(m0*width0*F,0.0);
313 double AA = m02-m2;
314 double BB = -m0*width;
315 EvtComplex amp = C/EvtComplex(AA,BB);
316 amplitude = abs(amp);
317 delta = atan2(imag(amp), real(amp));
318
319 double alpham2 = alpha*2.0;
320 F11 = amp*alpham2*q*H0*root2;
321 F21 = amp*alpham2*q*(Hp+Hm);
322 F31 = amp*alpham2*q*(Hp-Hm);
323}
324
325void EvtDToKpienu::NRS(double m, double q, double rS, double rS1, double a_delta, double b_delta, double mA, double m0, double width0, double& amplitude, double& delta, EvtComplex& F10)
326{
327 static const double tmp = (mK+mPi)*(mK+mPi);
328
329 double m2 = m*m;
330 double q2 = q*q;
331 double mA2 = mA*mA;
332 double pKPi = getPStar(mD, m, q);
333 double m_K0_1430 = m0;
334 double width_K0_1430 = width0;
335 double m2_K0_1430 = m_K0_1430*m_K0_1430;
336 double width = getWidth0(m, m_K0_1430, mPi, mK, width_K0_1430);
337
338 //calculate modul of the amplitude
339 double x,Pm;
340 if(m<m_K0_1430) {
341 x = sqrt(m2/tmp-1.0);
342 Pm = 1.0+rS1*x;
343 } else {
344 x = sqrt(m2_K0_1430/tmp-1.0);
345 Pm = 1.0+rS1*x;
346 Pm *= m_K0_1430*width_K0_1430/sqrt((m2_K0_1430-m2)*(m2_K0_1430-m2)+m2_K0_1430*width*width);
347 }
348
349 //calculate phase of the amplitude
350 double pStar = getPStar(m, mPi, mK);
351 double delta_bg = atan(2.*a_delta*pStar/(2.+a_delta*b_delta*pStar*pStar));
352 delta_bg = (delta_bg>0)? delta_bg: (delta_bg+Pi);
353
354 double delta_K0_1430 = atan(m_K0_1430*width/(m2_K0_1430-m2));
355 delta_K0_1430 = (delta_K0_1430>0)? delta_K0_1430: (delta_K0_1430+Pi);
356 delta = delta_bg + delta_K0_1430;
357
358 EvtComplex ci(cos(delta), sin(delta));
359 EvtComplex amp = ci*rS*Pm;
360 amplitude = rS*Pm;
361
362 F10 = amp*pKPi*mD/(1.-q2/mA2);
363}
364
365void EvtDToKpienu::ResonanceD(double m, double q, double mV, double mA, double TV_0, double T1_0, double T2_0, double m0, double width0, double rBW, double& amplitude, double& delta, EvtComplex& F12, EvtComplex& F22, EvtComplex& F32)
366{
367 double pKPi = getPStar(mD, m, q);
368 double mD2 = mD*mD;
369 double m2 = m*m;
370 double m02 = m0*m0;
371 double q2 = q*q;
372 double mV2 = mV*mV;
373 double mA2 = mA*mA;
374 double summDm = mD+m;
375 double TV = TV_0 /(1.0-q2/(mV2));
376 double T1 = T1_0/(1.0-q2/(mA2));
377 double T2 = T2_0/(1.0-q2/(mA2));
378
379 //construct amplitudes of (non)resonance
380 double F = getF2(m, m0, mPi, mK, rBW);
381 double width = getWidth2(m, m0, mPi, mK, width0, rBW);
382 EvtComplex C(m0*width0*F,0.0);
383 double AA = m02-m2;
384 double BB = -m0*width;
385
386 EvtComplex amp = C/EvtComplex(AA,BB);
387 amplitude = abs(amp);
388 delta = atan2(imag(amp), real(amp));
389
390 F12 = amp*mD*pKPi/3.*((mD2-m2-q2)*summDm*T1-mD2*pKPi*pKPi/summDm*T2);
391 F22 = amp*root2d3*mD*m*q*pKPi*summDm*T1;
392 F32 = amp*root2d3*2.*mD2*m*q*pKPi*pKPi/summDm*TV;
393}
394
395double EvtDToKpienu::getPStar(double m, double m1, double m2)
396{
397 double s = m*m;
398 double s1 = m1*m1;
399 double s2 = m2*m2;
400 double x = s + s1 - s2;
401 double t = 0.25*x*x/s - s1;
402 double p;
403 if (t>0.0) {
404 p = sqrt(t);
405 } else {
406 std::cout << " Hello, pstar is less than 0.0" << std::endl;
407 p = 0.04;
408 }
409 return p;
410}
411
412double EvtDToKpienu::getF1(double m, double m0, double m_c1, double m_c2, double rBW)
413{
414 double pStar = getPStar(m, m_c1, m_c2);
415 double pStar0 = getPStar(m0, m_c1, m_c2);
416 double rBW2 = rBW*rBW;
417 double pStar2 = pStar*pStar;
418 double pStar02 = pStar0*pStar0;
419 double B = 1./sqrt(1.+rBW2*pStar2);
420 double B0 = 1./sqrt(1.+rBW2*pStar02);
421 double F = pStar/pStar0*B/B0;
422 return F;
423}
424
425double EvtDToKpienu::getF2(double m, double m0, double m_c1, double m_c2, double rBW)
426{
427 double pStar = getPStar(m, m_c1, m_c2);
428 double pStar0 = getPStar(m0, m_c1, m_c2);
429 double rBW2 = rBW*rBW;
430 double pStar2 = pStar*pStar;
431 double pStar02 = pStar0*pStar0;
432 double B = 1./sqrt((rBW2*pStar2-3.)*(rBW2*pStar2-3.)+9.*rBW2*pStar2);
433 double B0 = 1./sqrt((rBW2*pStar02-3.)*(rBW2*pStar02-3.)+9.*rBW2*pStar02);
434 double F = pStar2/pStar02*B/B0;
435 return F;
436}
437
438double EvtDToKpienu::getWidth0(double m, double m0, double m_c1, double m_c2, double width0)
439{
440 double pStar = getPStar(m, m_c1, m_c2);
441 double pStar0 = getPStar(m0, m_c1, m_c2);
442 double width = width0*pStar/pStar0*m0/m;
443 return width;
444}
445
446double EvtDToKpienu::getWidth1(double m, double m0, double m_c1, double m_c2, double width0, double rBW)
447{
448 double pStar = getPStar(m, m_c1, m_c2);
449 double pStar0 = getPStar(m0, m_c1, m_c2);
450 double F = getF1(m, m0, m_c1, m_c2, rBW);
451 double width = width0*pStar/pStar0*m0/m*F*F;
452 return width;
453}
454
455double EvtDToKpienu::getWidth2(double m, double m0, double m_c1, double m_c2, double width0, double rBW)
456{
457 double pStar = getPStar(m, m_c1, m_c2);
458 double pStar0 = getPStar(m0, m_c1, m_c2);
459 double F = getF2(m, m0, m_c1, m_c2, rBW);
460 double width = width0*pStar/pStar0*m0/m*F*F;
461 return width;
462}
463
464EvtComplex EvtDToKpienu::getCoef(double rho, double phi)
465{
466 EvtComplex coef ( rho*cos(phi), rho*sin(phi) );
467 return coef;
468}
sin
double sin(const BesAngle a)
Definition: BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: BesAngle.h:213
delta
const double delta
Definition: CheckRung/check/gainVsDoca.cxx:32
x
Double_t x[10]
Definition: DataBase/tau_mode.c:57
I
const DifComplex I
alpha
const double alpha
Definition: FastVertexFit.cxx:4
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
C
***************************************************************************************Pseudo Class RRes *****************************************************************************************Parameters and physical constants **Maarten sept ************************************************************************DOUBLE PRECISION xsmu **************************************************************************PARTICLE DATA all others are from PDG *Only resonances with known widths into electron pairs are sept ************************************************************************C Declarations C
Definition: RRes.h:29
conj
Evt3Rank3C conj(const Evt3Rank3C &t2)
Definition: bak-BesEvtGen-00-04-08/src/EvtGen/EvtGenBase/Evt3Rank3C.cc:175
imag
double imag(const EvtComplex &c)
Definition: bak-BesEvtGen-00-04-08/src/EvtGen/EvtGenBase/EvtComplex.hh:246
abs2
double abs2(const EvtComplex &c)
Definition: bak-BesEvtGen-00-04-08/src/EvtGen/EvtGenBase/EvtComplex.hh:221
boostTo
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
Definition: bak-BesEvtGen-00-04-08/src/EvtGen/EvtGenBase/EvtDiracSpinor.cc:75
t
TTree * t
Definition: binning.cxx:23
EvtDecayBase::checkSpinDaughter
void checkSpinDaughter(int d1, EvtSpinType::spintype sp)
Definition: bak-BesEvtGen-00-04-08/src/EvtGen/EvtGenBase/EvtDecayBase.cc:533
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
pi
const float pi
Definition: vector3.h:133