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