49 <<
"EvtDToKpienu ==> Initialization !" << std::endl;
91 m_Pi = atan2( 0.0, -1.0 );
146 double m2, q2, cosV, cosL, chi;
147 KinVGen( K,
pi, e, nu, charm, m2, q2, cosV, cosL, chi );
148 double prob =
calPDF( m2, q2, cosV, cosL, chi );
155 const int charm,
double& m2,
double& q2,
156 double& cosV,
double& cosL,
double& chi )
const
161 m2 = vp4_KPi.
mass2();
165 boost.
set( vp4_KPi.
get( 0 ), -vp4_KPi.
get( 1 ), -vp4_KPi.
get( 2 ),
168 cosV = vp4_Kp.
dot( vp4_KPi ) / ( vp4_Kp.
d3mag() * vp4_KPi.
d3mag() );
170 boost.
set( vp4_W.
get( 0 ), -vp4_W.
get( 1 ), -vp4_W.
get( 2 ), -vp4_W.
get( 3 ) );
172 cosL = vp4_Lepp.
dot( vp4_W ) / ( vp4_Lepp.
d3mag() * vp4_W.
d3mag() );
179 double sinx = C.cross( V ).dot( D );
180 double cosx = C.dot( D );
181 chi = sinx > 0 ? acos( cosx ) : -acos( cosx );
187 const double cosL,
const double chi )
const
189 double m = sqrt( m2 );
190 double q = sqrt( q2 );
209 double amplitude_temp, delta_temp;
211 for (
int index = 0; index <
m_nAmps; index++ ) {
212 switch (
m_type[index] ) {
226 F11 = F11 + coef * f11;
227 F21 = F21 + coef * f21;
228 F31 = F31 + coef * f31;
237 F11 = F11 + coef * f11;
238 F21 = F21 + coef * f21;
239 F31 = F31 + coef * f31;
248 F12 = F12 + coef * f12;
249 F22 = F22 + coef * f22;
250 F32 = F32 + coef * f32;
255 <<
"No such form factor type!!!" << std::endl;
262 double I, I1, I2, I3, I4, I5, I6, I7, I8, I9;
264 double cosV2 = cosV * cosV;
265 double sinV = sqrt( 1.0 - cosV2 );
266 double sinV2 = sinV * sinV;
268 EvtComplex F1 = F10 + F11 * cosV + F12 * ( 1.5 * cosV2 - 0.5 );
272 I1 = 0.25 * (
abs2( F1 ) + 1.5 * sinV2 * (
abs2( F2 ) +
abs2( F3 ) ) );
273 I2 = -0.25 * (
abs2( F1 ) - 0.5 * sinV2 * (
abs2( F2 ) +
abs2( F3 ) ) );
274 I3 = -0.25 * (
abs2( F2 ) -
abs2( F3 ) ) * sinV2;
275 I4 =
real(
conj( F1 ) * F2 ) * sinV * 0.5;
276 I5 =
real(
conj( F1 ) * F3 ) * sinV;
277 I6 =
real(
conj( F2 ) * F3 ) * sinV2;
278 I7 =
imag(
conj( F2 ) * F1 ) * sinV;
279 I8 =
imag(
conj( F3 ) * F1 ) * sinV * 0.5;
280 I9 =
imag(
conj( F3 ) * F2 ) * sinV2 * ( -0.5 );
282 double coschi = cos( chi );
283 double sinchi = sin( chi );
284 double sin2chi = 2.0 * sinchi * coschi;
285 double cos2chi = 1.0 - 2.0 * sinchi * sinchi;
287 double sinL = sqrt( 1. - cosL * cosL );
288 double sinL2 = sinL * sinL;
289 double sin2L = 2.0 * sinL * cosL;
290 double cos2L = 1.0 - 2.0 * sinL2;
292 I = I1 + I2 * cos2L + I3 * sinL2 * cos2chi + I4 * sin2L * coschi +
293 I5 * sinL * coschi + I6 * cosL + I7 * sinL * sinchi +
294 I8 * sin2L * sinchi + I9 * sinL2 * sin2chi;
299 const double mA,
const double V_0,
300 const double A1_0,
const double A2_0,
301 const double m0,
const double width0,
302 const double rBW,
double& amplitude,
309 double m02 = m0 * m0;
311 double mV2 = mV * mV;
312 double mA2 = mA * mA;
313 double summDm =
m_mD + m;
314 double V = V_0 / ( 1.0 - q2 / ( mV2 ) );
315 double A1 = A1_0 / ( 1.0 - q2 / ( mA2 ) );
316 double A2 = A2_0 / ( 1.0 - q2 / ( mA2 ) );
317 double A = summDm * A1;
318 double B = 2.0 *
m_mD * pKPi / summDm * V;
321 double H0 = 0.5 / ( m * q ) *
322 ( ( mD2 - m2 - q2 ) * summDm * A1 -
323 4.0 * ( mD2 * pKPi * pKPi ) / summDm * A2 );
328 double B_Kstar = 2. / 3.;
330 double alpha = sqrt( 3. *
m_Pi * B_Kstar / ( pStar0 * width0 ) );
337 double AA = m02 - m2;
338 double BB = -m0 * width;
340 amplitude =
abs( amp );
341 delta = atan2(
imag( amp ),
real( amp ) );
343 double alpham2 = alpha * 2.0;
344 F11 = amp * alpham2 * q * H0 *
m_root2;
345 F21 = amp * alpham2 * q * ( Hp + Hm );
346 F31 = amp * alpham2 * q * ( Hp - Hm );
350 const double rS1,
const double a_delta,
351 const double b_delta,
const double mA,
const double m0,
352 const double width0,
double& amplitude,
double& delta,
359 double mA2 = mA * mA;
361 double m_K0_1430 = m0;
362 double width_K0_1430 = width0;
363 double m2_K0_1430 = m_K0_1430 * m_K0_1430;
368 if ( m < m_K0_1430 ) {
369 x = sqrt( m2 / tmp - 1.0 );
372 x = sqrt( m2_K0_1430 / tmp - 1.0 );
374 Pm *= m_K0_1430 * width_K0_1430 /
375 sqrt( ( m2_K0_1430 - m2 ) * ( m2_K0_1430 - m2 ) +
376 m2_K0_1430 * width * width );
381 double delta_bg = atan( 2. * a_delta * pStar /
382 ( 2. + a_delta * b_delta * pStar * pStar ) );
383 delta_bg = ( delta_bg > 0 ) ? delta_bg : ( delta_bg +
m_Pi );
385 double delta_K0_1430 = atan( m_K0_1430 * width / ( m2_K0_1430 - m2 ) );
386 delta_K0_1430 = ( delta_K0_1430 > 0 ) ? delta_K0_1430
387 : ( delta_K0_1430 +
m_Pi );
388 delta = delta_bg + delta_K0_1430;
394 F10 = amp * pKPi *
m_mD / ( 1. - q2 / mA2 );
398 const double mA,
const double TV_0,
399 const double T1_0,
const double T2_0,
400 const double m0,
const double width0,
401 const double rBW,
double& amplitude,
408 double m02 = m0 * m0;
410 double mV2 = mV * mV;
411 double mA2 = mA * mA;
412 double summDm =
m_mD + m;
413 double TV = TV_0 / ( 1.0 - q2 / ( mV2 ) );
414 double T1 = T1_0 / ( 1.0 - q2 / ( mA2 ) );
415 double T2 = T2_0 / ( 1.0 - q2 / ( mA2 ) );
421 double AA = m02 - m2;
422 double BB = -m0 * width;
425 amplitude =
abs( amp );
426 delta = atan2(
imag( amp ),
real( amp ) );
428 F12 = amp *
m_mD * pKPi / 3. *
429 ( ( mD2 - m2 - q2 ) * summDm * T1 - mD2 * pKPi * pKPi / summDm * T2 );
431 F32 = amp *
m_root2d3 * 2. * mD2 * m * q * pKPi * pKPi / summDm * TV;
435 const double m2 )
const
440 double x = s + s1 - s2;
441 double t = 0.25 * x * x / s - s1;
447 <<
" Hello, pstar is less than 0.0" << std::endl;
454 const double m_c2,
const double rBW )
const
456 double pStar =
getPStar( m, m_c1, m_c2 );
457 double pStar0 =
getPStar( m0, m_c1, m_c2 );
458 double rBW2 = rBW * rBW;
459 double pStar2 = pStar * pStar;
460 double pStar02 = pStar0 * pStar0;
461 double B = 1. / sqrt( 1. + rBW2 * pStar2 );
462 double B0 = 1. / sqrt( 1. + rBW2 * pStar02 );
463 double F = pStar / pStar0 * B / B0;
468 const double m_c2,
const double rBW )
const
470 double pStar =
getPStar( m, m_c1, m_c2 );
471 double pStar0 =
getPStar( m0, m_c1, m_c2 );
472 double rBW2 = rBW * rBW;
473 double pStar2 = pStar * pStar;
474 double pStar02 = pStar0 * pStar0;
475 double B = 1. / sqrt( ( rBW2 * pStar2 - 3. ) * ( rBW2 * pStar2 - 3. ) +
476 9. * rBW2 * pStar2 );
477 double B0 = 1. / sqrt( ( rBW2 * pStar02 - 3. ) * ( rBW2 * pStar02 - 3. ) +
478 9. * rBW2 * pStar02 );
479 double F = pStar2 / pStar02 * B / B0;
484 const double m_c1,
const double m_c2,
485 const double width0 )
const
487 double pStar =
getPStar( m, m_c1, m_c2 );
488 double pStar0 =
getPStar( m0, m_c1, m_c2 );
489 double width = width0 * pStar / pStar0 * m0 / m;
494 const double m_c1,
const double m_c2,
495 const double width0,
const double rBW )
const
497 double pStar =
getPStar( m, m_c1, m_c2 );
498 double pStar0 =
getPStar( m0, m_c1, m_c2 );
499 double F =
getF1( m, m0, m_c1, m_c2, rBW );
500 double width = width0 * pStar / pStar0 * m0 / m * F * F;
505 const double m_c1,
const double m_c2,
506 const double width0,
const double rBW )
const
508 double pStar =
getPStar( m, m_c1, m_c2 );
509 double pStar0 =
getPStar( m0, m_c1, m_c2 );
510 double F =
getF2( m, m0, m_c1, m_c2, rBW );
511 double width = width0 * pStar / pStar0 * m0 / m * F * F;
517 EvtComplex coef( rho * cos( phi ), rho * sin( phi ) );
EvtComplex conj(const EvtComplex &c)
double imag(const EvtComplex &c)
double real(const EvtComplex &c)
double abs2(const EvtComplex &c)
double abs(const EvtComplex &c)
EvtRaritaSchwinger boostTo(const EvtRaritaSchwinger &rs, const EvtVector4R p4)
std::ostream & EvtGenReport(EvtGenSeverity severity, const char *facility=nullptr)
double getWidth0(const double m, const double m0, const double m_c1, const double m_c2, const double width0) const
double getPStar(const double m, const double m1, const double m2) const
void initProbMax() override
void KinVGen(const EvtVector4R &vp4_K, const EvtVector4R &vp4_Pi, const EvtVector4R &vp4_Lep, const EvtVector4R &vp4_Nu, const int charm, double &m2, double &q2, double &cosV, double &cosL, double &chi) const
void ResonanceD(const double m, const double q, const double mV, const double mA, const double TV_0, const double T1_0, const double T2_0, const double m0, const double width0, const double rBW, double &litude, double &delta, EvtComplex &F12, EvtComplex &F22, EvtComplex &F32) const
double getWidth2(const double m, const double m0, const double m_c1, const double m_c2, const double width0, const double rBW) const
EvtDecayBase * clone() const override
double getF2(const double m, const double m0, const double m_c1, const double m_c2, const double rBW) const
double getWidth1(const double m, const double m0, const double m_c1, const double m_c2, const double width0, const double rBW) const
void ResonanceP(const double m, const double q, const double mV, const double mA, const double V_0, const double A1_0, const double A2_0, const double m0, const double width0, const double rBW, double &litude, double &delta, EvtComplex &F11, EvtComplex &F21, EvtComplex &F31) const
void decay(EvtParticle *p) override
void NRS(const double m, const double q, const double rS, const double rS1, const double a_delta, const double b_delta, const double mA, const double m0, const double width0, double &litude, double &delta, EvtComplex &F10) const
std::array< int, 5 > m_type
double getF1(const double m, const double m0, const double m_c1, const double m_c2, const double rBW) const
double calPDF(const double m2, const double q2, const double cosV, const double cosL, const double chi) const
EvtComplex getCoef(const double rho, const double phi) const
std::string getName() const override
void checkSpinDaughter(int d1, EvtSpinType::spintype sp)
void checkSpinParent(EvtSpinType::spintype sp)
void setProbMax(double prbmx)
void checkNDaug(int d1, int d2=-1)
void checkNArg(int a1, int a2=-1, int a3=-1, int a4=-1)
const EvtId * getDaugs() const
void setProb(double prob)
static int getStdHep(EvtId id)
double initializePhaseSpace(size_t numdaughter, const EvtId *daughters, bool forceResetMasses=false, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)
const EvtVector4R & getP4() const
EvtParticle * getDaug(const int i)
double dot(const EvtVector4R &v2) const
EvtVector4R cross(const EvtVector4R &v2) const
void set(int i, double d)