47 tan2thetaC = ( 0.22650 * 0.22650 ) /
48 ( 1. - ( 0.22650 * 0.22650 ) );
101 mass_R[10] = 1.41400;
102 width_R[10] = 0.2320;
106 mass_R[11] = 1.42500;
107 width_R[11] = 0.2700;
111 mass_R[12] = 1.42500;
112 width_R[12] = 0.2700;
118 EvtComplex( 0.255303 *
cos( 47.8861 * pi180inv ), 0.255303 *
sin( 47.8861 * pi180inv ) );
120 EvtComplex( 13.4446 *
cos( -5.11127 * pi180inv ), 13.4446 *
sin( -5.11127 * pi180inv ) );
122 EvtComplex( 38.8496 *
cos( -30.06 * pi180inv ), 38.8496 *
sin( -30.06 * pi180inv ) );
124 EvtComplex( 13.1086 *
cos( -81.4148 * pi180inv ), 13.1086 *
sin( -81.4148 * pi180inv ) );
128 EvtComplex( 5.08049 *
cos( -182.312 * pi180inv ), 5.08049 *
sin( -182.312 * pi180inv ) );
130 EvtComplex( 17.2388 *
cos( -219.209 * pi180inv ), 17.2388 *
sin( -219.209 * pi180inv ) );
132 EvtComplex( 19.0145 *
cos( -76.9884 * pi180inv ), 19.0145 *
sin( -76.9884 * pi180inv ) );
134 EvtComplex( 11.9875 *
cos( -190.502 * pi180inv ), 11.9875 *
sin( -190.502 * pi180inv ) );
174 deltad[1] = 194.7 * pi180inv;
175 deltad[2] = 196.0 * pi180inv;
176 deltad[3] = 167.0 * pi180inv;
184 double ProbMax = 12500.0;
185 if ( tagmode == 1 ) ProbMax = 12500.0;
186 else if ( tagmode == 2 ) ProbMax = 12000.0;
187 else if ( tagmode == 3 ) ProbMax = 12100.0;
188 else ProbMax = 12000.0;
210 for (
int i = 0; i < _nd; i++ )
215 double prob = AmplitudeSquare();
220double EvtD0ToKSpipi::AmplitudeSquare() {
226 double total_prob, Interference = 0.;
229 EvtComplex DK2piRes0 = Resonance2( pD, pip, pim, ar[0], phir[0], width_R[0], mass_R[0],
231 EvtComplex DK2piRes1 = Resonance2( pD, pip, pim, ar[1], phir[1], width_R[1], mass_R[1],
233 EvtComplex DK2piRes2 = Resonance2( pD, pip, pim, ar[2], phir[2], width_R[2], mass_R[2],
235 EvtComplex DK2piRes3 = Resonance2( pD, pip, pim, ar[3], phir[3], width_R[3], mass_R[3],
237 EvtComplex DK2piRes4 = Resonance2( pD, k0l, pim, ar[4], phir[4], width_R[4], mass_R[4],
239 EvtComplex DK2piRes5 = Resonance2( pD, k0l, pim, ar[5], phir[5], width_R[5], mass_R[5],
241 EvtComplex DK2piRes6 = Resonance2( pD, k0l, pim, ar[6], phir[6], width_R[6], mass_R[6],
243 EvtComplex DK2piRes7 = Resonance2( pD, k0l, pim, ar[7], phir[7], width_R[7], mass_R[7],
245 EvtComplex DK2piRes8 = Resonance2( pD, k0l, pip, ar[8], phir[8], width_R[8], mass_R[8],
247 EvtComplex DK2piRes9 = Resonance2( pD, k0l, pip, ar[9], phir[9], width_R[9], mass_R[9],
249 EvtComplex DK2piRes10 = Resonance2( pD, k0l, pip, ar[10], phir[10], width_R[10], mass_R[10],
253 EvtComplex pipi_s_wave = K_matrix( pip, pim );
254 if ( pipi_s_wave ==
EvtComplex( 9999., 9999. ) )
return 1e-20;
258 amplitude_LASS( k0l, pip, pim,
"k0spim", ar[11], phir[11] * pi180inv );
260 amplitude_LASS( k0l, pip, pim,
"k0spip", ar[12], phir[12] * pi180inv );
263 EvtComplex TOT_PFT_AMP = DK2piRes0 + DK2piRes1 + DK2piRes2 + DK2piRes3 + DK2piRes4 +
264 DK2piRes5 + DK2piRes6 + DK2piRes7 + DK2piRes8 + DK2piRes9 +
265 DK2piRes10 + pipi_s_wave + kpi_s_wave + dcs_kpi_s_wave;
267 if ( tagmode == 1 || tagmode == 2 || tagmode == 3 )
270 EvtComplex DK2piRes4_dcs = Resonance2( pD, k0l, pim, ar[8], phir[8], width_R[4], mass_R[4],
272 EvtComplex DK2piRes5_dcs = Resonance2( pD, k0l, pim, ar[9], phir[9], width_R[5], mass_R[5],
274 EvtComplex DK2piRes6_dcs = Resonance2( pD, k0l, pip, ar[6], phir[6], width_R[6], mass_R[6],
276 EvtComplex DK2piRes7_dcs = Resonance2( pD, k0l, pim, ar[10], phir[10], width_R[7],
277 mass_R[7], spin_R[7] );
278 EvtComplex DK2piRes8_dcs = Resonance2( pD, k0l, pip, ar[4], phir[4], width_R[8], mass_R[8],
280 EvtComplex DK2piRes9_dcs = Resonance2( pD, k0l, pip, ar[5], phir[5], width_R[9], mass_R[9],
282 EvtComplex DK2piRes10_dcs = Resonance2( pD, k0l, pip, ar[7], phir[7], width_R[10],
283 mass_R[10], spin_R[10] );
286 amplitude_LASS( k0l, pip, pim,
"k0spip", ar[11], phir[11] * pi180inv );
288 amplitude_LASS( k0l, pip, pim,
"k0spim", ar[12], phir[12] * pi180inv );
291 EvtComplex TOT_DCS_AMP = DK2piRes0 + DK2piRes1 + DK2piRes2 + DK2piRes3 + DK2piRes4_dcs +
292 DK2piRes5_dcs + DK2piRes6_dcs + DK2piRes7_dcs + DK2piRes8_dcs +
293 DK2piRes9_dcs + DK2piRes10_dcs + pipi_s_wave + kpi_s_wave_dcs +
299 TOT_DCS_AMP *
conj( TOT_PFT_AMP ) +
301 TOT_PFT_AMP *
conj( TOT_DCS_AMP ) );
304 total_prob =
abs2( TOT_PFT_AMP ) + rd[tagmode] * rd[tagmode] *
abs2( TOT_DCS_AMP ) -
305 rd[tagmode] * Rf[tagmode] * Interference;
307 else { total_prob =
abs2( TOT_PFT_AMP ); }
309 return ( total_prob <= 0 ) ? 1e-20 : total_prob;
314 double gamma,
double bwm,
int spin ) {
317 EvtVector4R p4_d3 = p4_p - p4_d1 - p4_d2;
319 double mAB = ( p4_d1 + p4_d2 ).
mass();
320 double mBC = ( p4_d2 + p4_d3 ).
mass();
321 double mAC = ( p4_d1 + p4_d3 ).
mass();
322 double mA = p4_d1.
mass();
323 double mB = p4_d2.
mass();
324 double mD = p4_p.
mass();
325 double mC = p4_d3.
mass();
328 double gammaR = gamma;
330 sqrt( ( ( ( mAB * mAB - mA * mA - mB * mB ) * ( mAB * mAB - mA * mA - mB * mB ) / 4.0 ) -
331 mA * mA * mB * mB ) /
334 sqrt( ( ( ( mR * mR - mA * mA - mB * mB ) * ( mR * mR - mA * mA - mB * mB ) / 4.0 ) -
335 mA * mA * mB * mB ) /
338 double pD = ( ( ( mD * mD - mR * mR - mC * mC ) * ( mD * mD - mR * mR - mC * mC ) / 4.0 ) -
339 mR * mR * mC * mC ) /
341 if ( pD > 0 ) { pD = sqrt( pD ); }
344 sqrt( ( ( ( mD * mD - mAB * mAB - mC * mC ) * ( mD * mD - mAB * mAB - mC * mC ) / 4.0 ) -
345 mAB * mAB * mC * mC ) /
358 fR = sqrt( 1.0 + 1.5 * 1.5 * pR * pR ) / sqrt( 1.0 + 1.5 * 1.5 * pAB * pAB );
359 fD = sqrt( 1.0 + 5.0 * 5.0 * pD * pD ) / sqrt( 1.0 + 5.0 * 5.0 * pDAB * pDAB );
363 fR = sqrt( ( 9 + 3 * pow( ( 1.5 * pR ), 2 ) + pow( ( 1.5 * pR ), 4 ) ) /
364 ( 9 + 3 * pow( ( 1.5 * pAB ), 2 ) + pow( ( 1.5 * pAB ), 4 ) ) );
365 fD = sqrt( ( 9 + 3 * pow( ( 5.0 * pD ), 2 ) + pow( ( 5.0 * pD ), 4 ) ) /
366 ( 9 + 3 * pow( ( 5.0 * pDAB ), 2 ) + pow( ( 5.0 * pDAB ), 4 ) ) );
369 default:
report(
INFO,
"EvtGen" ) <<
"Incorrect spin in EvtD0ToKSpipi::EvtResonance2.cc\n";
372 double gammaAB = gammaR * pow( pAB / pR, power ) * ( mR / mAB ) * fR * fR;
376 ampl = mag * EvtComplex(
cos( theta * pi180inv ),
sin( theta * pi180inv ) ) * fR * fD /
377 ( mR * mR - mAB * mAB - EvtComplex( 0.0, mR * gammaAB ) );
380 ampl = mag * EvtComplex(
cos( theta * pi180inv ),
sin( theta * pi180inv ) ) *
382 ( mAC * mAC - mBC * mBC +
383 ( ( mD * mD - mC * mC ) * ( mB * mB - mA * mA ) / ( mAB * mAB ) ) ) /
384 ( mR * mR - mAB * mAB - EvtComplex( 0.0, mR * gammaAB ) ) );
387 ampl = mag * EvtComplex(
cos( theta * pi180inv ),
sin( theta * pi180inv ) ) *
388 ( fR * fD / ( mR * mR - mAB * mAB - EvtComplex( 0.0, mR * gammaAB ) ) ) *
389 ( pow( ( mBC * mBC - mAC * mAC +
390 ( mD * mD - mC * mC ) * ( mA * mA - mB * mB ) / ( mAB * mAB ) ),
393 ( mAB * mAB - 2 * mD * mD - 2 * mC * mC +
394 pow( ( mD * mD - mC * mC ) / mAB, 2 ) ) *
395 ( mAB * mAB - 2 * mA * mA - 2 * mB * mB +
396 pow( ( mA * mA - mB * mB ) / mAB, 2 ) ) );
398 default:
report(
INFO,
"EvtGen" ) <<
"Incorrect spin in EvtD0ToKSpipi::Resonance2.cc\n";
405 const double mD0 = 1.86483;
406 const double mKl = 0.49761;
407 const double mPi = 0.13957;
410 double mAB = ( p_pip + p_pim ).
mass();
411 double s = mAB * mAB;
413 EvtComplex n11, n12, n13, n14, n15, n21, n22, n23, n24, n25, n31, n32, n33, n34, n35, n41,
414 n42, n43, n44, n45, n51, n52, n53, n54, n55;
415 double rho1sq, rho2sq, rho4sq, rho5sq;
416 EvtComplex rho1, rho2, rho3, rho4, rho5;
418 EvtComplex pole, SVT, Adler;
423 const double mpi = 0.13957;
424 const double mK = 0.493677;
425 const double meta = 0.54775;
426 const double metap = 0.95778;
430 for (
int k = 0; k < 5; k++ )
432 for (
int l = 0; l < 5; l++ )
434 i[k][l] = EvtComplex( 0., 0. );
435 K[k][l] = EvtComplex( 0., 0. );
442 double s_scatt = -3.92637;
460 rho1sq = ( 1.0 - ( pow( (
mpi +
mpi ), 2 ) /
s ) );
461 if ( rho1sq >= 0. ) rho1 = EvtComplex( sqrt( rho1sq ), 0. );
462 else rho1 = EvtComplex( 0., sqrt( -rho1sq ) );
466 rho2sq = ( 1.0 - ( pow( ( mK + mK ), 2 ) /
s ) );
467 if ( rho2sq >= 0. ) rho2 = EvtComplex( sqrt( rho2sq ), 0. );
468 else rho2 = EvtComplex( 0., sqrt( -rho2sq ) );
472 rho3 = EvtComplex( 0., 0. );
475 double real = 1.2274 + 0.00370909 / (
s *
s ) - ( 0.111203 ) / (
s)-6.39017 *
s +
476 16.8358 *
s *
s - 21.8845 *
s *
s *
s + 11.3153 *
s *
s *
s *
s;
477 double cont32 = sqrt( 1.0 - ( 16.0 *
mpi *
mpi ) );
478 rho3 = EvtComplex( cont32 *
real, 0. );
480 else rho3 = EvtComplex( sqrt( 1.0 - ( 16.0 *
mpi *
mpi /
s ) ), 0. );
484 rho4sq = ( 1.0 - ( pow( (
meta +
meta ), 2 ) /
s ) );
485 if ( rho4sq >= 0. ) rho4 = EvtComplex( sqrt( rho4sq ), 0. );
486 else rho4 = EvtComplex( 0., sqrt( -rho4sq ) );
490 rho5sq = ( 1.0 - ( pow( (
meta + metap ), 2 ) /
s ) );
491 if ( rho5sq >= 0. ) rho5 = EvtComplex( sqrt( rho5sq ), 0. );
492 else rho5 = EvtComplex( 0., sqrt( -rho5sq ) );
496 for (
int k = 0; k < 5; k++ )
498 for (
int l = 0; l < 5; l++ )
500 for (
int pole_index = 0; pole_index < 5; pole_index++ )
502 double A = g[pole_index][k] * g[pole_index][l];
503 double B = ma[pole_index] * ma[pole_index] -
s;
504 K[k][l] = K[k][l] + EvtComplex( A / B, 0. );
510 for (
int k = 0; k < 5; k++ )
512 for (
int l = 0; l < 5; l++ )
514 double C =
f[k][l] * ( 1.0 - s_scatt );
515 double D = (
s - s_scatt );
516 K[k][l] = K[k][l] + EvtComplex(
C / D, 0. );
521 for (
int k = 0; k < 5; k++ )
523 for (
int l = 0; l < 5; l++ )
525 double E = (
s - ( sa *
mpi *
mpi * 0.5 ) ) * ( 1.0 - sa_0 );
526 double F = (
s - sa_0 );
527 K[k][l] = K[k][l] * EvtComplex( E / F, 0. );
532 n11 = EvtComplex( 1., 0. ) - EvtComplex( 0., 1. ) * K[0][0] * rho[0];
533 n12 = EvtComplex( 0., 0. ) - EvtComplex( 0., 1. ) * K[0][1] * rho[1];
534 n13 = EvtComplex( 0., 0. ) - EvtComplex( 0., 1. ) * K[0][2] * rho[2];
535 n14 = EvtComplex( 0., 0. ) - EvtComplex( 0., 1. ) * K[0][3] * rho[3];
536 n15 = EvtComplex( 0., 0. ) - EvtComplex( 0., 1. ) * K[0][4] * rho[4];
538 n21 = EvtComplex( 0., 0. ) - EvtComplex( 0., 1. ) * K[1][0] * rho[0];
539 n22 = EvtComplex( 1., 0. ) - EvtComplex( 0., 1. ) * K[1][1] * rho[1];
540 n23 = EvtComplex( 0., 0. ) - EvtComplex( 0., 1. ) * K[1][2] * rho[2];
541 n24 = EvtComplex( 0., 0. ) - EvtComplex( 0., 1. ) * K[1][3] * rho[3];
542 n25 = EvtComplex( 0., 0. ) - EvtComplex( 0., 1. ) * K[1][4] * rho[4];
544 n31 = EvtComplex( 0., 0. ) - EvtComplex( 0., 1. ) * K[2][0] * rho[0];
545 n32 = EvtComplex( 0., 0. ) - EvtComplex( 0., 1. ) * K[2][1] * rho[1];
546 n33 = EvtComplex( 1., 0. ) - EvtComplex( 0., 1. ) * K[2][2] * rho[2];
547 n34 = EvtComplex( 0., 0. ) - EvtComplex( 0., 1. ) * K[2][3] * rho[3];
548 n35 = EvtComplex( 0., 0. ) - EvtComplex( 0., 1. ) * K[2][4] * rho[4];
550 n41 = EvtComplex( 0., 0. ) - EvtComplex( 0., 1. ) * K[3][0] * rho[0];
551 n42 = EvtComplex( 0., 0. ) - EvtComplex( 0., 1. ) * K[3][1] * rho[1];
552 n43 = EvtComplex( 0., 0. ) - EvtComplex( 0., 1. ) * K[3][2] * rho[2];
553 n44 = EvtComplex( 1., 0. ) - EvtComplex( 0., 1. ) * K[3][3] * rho[3];
554 n45 = EvtComplex( 0., 0. ) - EvtComplex( 0., 1. ) * K[3][4] * rho[4];
556 n51 = EvtComplex( 0., 0. ) - EvtComplex( 0., 1. ) * K[4][0] * rho[0];
557 n52 = EvtComplex( 0., 0. ) - EvtComplex( 0., 1. ) * K[4][1] * rho[1];
558 n53 = EvtComplex( 0., 0. ) - EvtComplex( 0., 1. ) * K[4][2] * rho[2];
559 n54 = EvtComplex( 0., 0. ) - EvtComplex( 0., 1. ) * K[4][3] * rho[3];
560 n55 = EvtComplex( 1., 0. ) - EvtComplex( 0., 1. ) * K[4][4] * rho[4];
564 det = ( n15 * n24 * n33 * n42 * n51 - n14 * n25 * n33 * n42 * n51 -
565 n15 * n23 * n34 * n42 * n51 + n13 * n25 * n34 * n42 * n51 +
566 n14 * n23 * n35 * n42 * n51 - n13 * n24 * n35 * n42 * n51 -
567 n15 * n24 * n32 * n43 * n51 + n14 * n25 * n32 * n43 * n51 +
568 n15 * n22 * n34 * n43 * n51 - n12 * n25 * n34 * n43 * n51 -
569 n14 * n22 * n35 * n43 * n51 + n12 * n24 * n35 * n43 * n51 +
570 n15 * n23 * n32 * n44 * n51 - n13 * n25 * n32 * n44 * n51 -
571 n15 * n22 * n33 * n44 * n51 + n12 * n25 * n33 * n44 * n51 +
572 n13 * n22 * n35 * n44 * n51 - n12 * n23 * n35 * n44 * n51 -
573 n14 * n23 * n32 * n45 * n51 + n13 * n24 * n32 * n45 * n51 +
574 n14 * n22 * n33 * n45 * n51 - n12 * n24 * n33 * n45 * n51 -
575 n13 * n22 * n34 * n45 * n51 + n12 * n23 * n34 * n45 * n51 -
576 n15 * n24 * n33 * n41 * n52 + n14 * n25 * n33 * n41 * n52 +
577 n15 * n23 * n34 * n41 * n52 - n13 * n25 * n34 * n41 * n52 -
578 n14 * n23 * n35 * n41 * n52 + n13 * n24 * n35 * n41 * n52 +
579 n15 * n24 * n31 * n43 * n52 - n14 * n25 * n31 * n43 * n52 -
580 n15 * n21 * n34 * n43 * n52 + n11 * n25 * n34 * n43 * n52 +
581 n14 * n21 * n35 * n43 * n52 - n11 * n24 * n35 * n43 * n52 -
582 n15 * n23 * n31 * n44 * n52 + n13 * n25 * n31 * n44 * n52 +
583 n15 * n21 * n33 * n44 * n52 - n11 * n25 * n33 * n44 * n52 -
584 n13 * n21 * n35 * n44 * n52 + n11 * n23 * n35 * n44 * n52 +
585 n14 * n23 * n31 * n45 * n52 - n13 * n24 * n31 * n45 * n52 -
586 n14 * n21 * n33 * n45 * n52 + n11 * n24 * n33 * n45 * n52 +
587 n13 * n21 * n34 * n45 * n52 - n11 * n23 * n34 * n45 * n52 +
588 n15 * n24 * n32 * n41 * n53 - n14 * n25 * n32 * n41 * n53 -
589 n15 * n22 * n34 * n41 * n53 + n12 * n25 * n34 * n41 * n53 +
590 n14 * n22 * n35 * n41 * n53 - n12 * n24 * n35 * n41 * n53 -
591 n15 * n24 * n31 * n42 * n53 + n14 * n25 * n31 * n42 * n53 +
592 n15 * n21 * n34 * n42 * n53 - n11 * n25 * n34 * n42 * n53 -
593 n14 * n21 * n35 * n42 * n53 + n11 * n24 * n35 * n42 * n53 +
594 n15 * n22 * n31 * n44 * n53 - n12 * n25 * n31 * n44 * n53 -
595 n15 * n21 * n32 * n44 * n53 + n11 * n25 * n32 * n44 * n53 +
596 n12 * n21 * n35 * n44 * n53 - n11 * n22 * n35 * n44 * n53 -
597 n14 * n22 * n31 * n45 * n53 + n12 * n24 * n31 * n45 * n53 +
598 n14 * n21 * n32 * n45 * n53 - n11 * n24 * n32 * n45 * n53 -
599 n12 * n21 * n34 * n45 * n53 + n11 * n22 * n34 * n45 * n53 -
600 n15 * n23 * n32 * n41 * n54 + n13 * n25 * n32 * n41 * n54 +
601 n15 * n22 * n33 * n41 * n54 - n12 * n25 * n33 * n41 * n54 -
602 n13 * n22 * n35 * n41 * n54 + n12 * n23 * n35 * n41 * n54 +
603 n15 * n23 * n31 * n42 * n54 - n13 * n25 * n31 * n42 * n54 -
604 n15 * n21 * n33 * n42 * n54 + n11 * n25 * n33 * n42 * n54 +
605 n13 * n21 * n35 * n42 * n54 - n11 * n23 * n35 * n42 * n54 -
606 n15 * n22 * n31 * n43 * n54 + n12 * n25 * n31 * n43 * n54 +
607 n15 * n21 * n32 * n43 * n54 - n11 * n25 * n32 * n43 * n54 -
608 n12 * n21 * n35 * n43 * n54 + n11 * n22 * n35 * n43 * n54 +
609 n13 * n22 * n31 * n45 * n54 - n12 * n23 * n31 * n45 * n54 -
610 n13 * n21 * n32 * n45 * n54 + n11 * n23 * n32 * n45 * n54 +
611 n12 * n21 * n33 * n45 * n54 - n11 * n22 * n33 * n45 * n54 +
612 n14 * n23 * n32 * n41 * n55 - n13 * n24 * n32 * n41 * n55 -
613 n14 * n22 * n33 * n41 * n55 + n12 * n24 * n33 * n41 * n55 +
614 n13 * n22 * n34 * n41 * n55 - n12 * n23 * n34 * n41 * n55 -
615 n14 * n23 * n31 * n42 * n55 + n13 * n24 * n31 * n42 * n55 +
616 n14 * n21 * n33 * n42 * n55 - n11 * n24 * n33 * n42 * n55 -
617 n13 * n21 * n34 * n42 * n55 + n11 * n23 * n34 * n42 * n55 +
618 n14 * n22 * n31 * n43 * n55 - n12 * n24 * n31 * n43 * n55 -
619 n14 * n21 * n32 * n43 * n55 + n11 * n24 * n32 * n43 * n55 +
620 n12 * n21 * n34 * n43 * n55 - n11 * n22 * n34 * n43 * n55 -
621 n13 * n22 * n31 * n44 * n55 + n12 * n23 * n31 * n44 * n55 +
622 n13 * n21 * n32 * n44 * n55 - n11 * n23 * n32 * n44 * n55 -
623 n12 * n21 * n33 * n44 * n55 + n11 * n22 * n33 * n44 * n55 );
625 if ( det == EvtComplex( 0., 0. ) ) reject =
true;
628 i[0][0] = ( n25 * n34 * n43 * n52 - n24 * n35 * n43 * n52 - n25 * n33 * n44 * n52 +
629 n23 * n35 * n44 * n52 + n24 * n33 * n45 * n52 - n23 * n34 * n45 * n52 -
630 n25 * n34 * n42 * n53 + n24 * n35 * n42 * n53 + n25 * n32 * n44 * n53 -
631 n22 * n35 * n44 * n53 - n24 * n32 * n45 * n53 + n22 * n34 * n45 * n53 +
632 n25 * n33 * n42 * n54 - n23 * n35 * n42 * n54 - n25 * n32 * n43 * n54 +
633 n22 * n35 * n43 * n54 + n23 * n32 * n45 * n54 - n22 * n33 * n45 * n54 -
634 n24 * n33 * n42 * n55 + n23 * n34 * n42 * n55 + n24 * n32 * n43 * n55 -
635 n22 * n34 * n43 * n55 - n23 * n32 * n44 * n55 + n22 * n33 * n44 * n55 ) /
638 i[0][1] = ( -n15 * n34 * n43 * n52 + n14 * n35 * n43 * n52 + n15 * n33 * n44 * n52 -
639 n13 * n35 * n44 * n52 - n14 * n33 * n45 * n52 + n13 * n34 * n45 * n52 +
640 n15 * n34 * n42 * n53 - n14 * n35 * n42 * n53 - n15 * n32 * n44 * n53 +
641 n12 * n35 * n44 * n53 + n14 * n32 * n45 * n53 - n12 * n34 * n45 * n53 -
642 n15 * n33 * n42 * n54 + n13 * n35 * n42 * n54 + n15 * n32 * n43 * n54 -
643 n12 * n35 * n43 * n54 - n13 * n32 * n45 * n54 + n12 * n33 * n45 * n54 +
644 n14 * n33 * n42 * n55 - n13 * n34 * n42 * n55 - n14 * n32 * n43 * n55 +
645 n12 * n34 * n43 * n55 + n13 * n32 * n44 * n55 - n12 * n33 * n44 * n55 ) /
648 i[0][2] = ( n15 * n24 * n43 * n52 - n14 * n25 * n43 * n52 - n15 * n23 * n44 * n52 +
649 n13 * n25 * n44 * n52 + n14 * n23 * n45 * n52 - n13 * n24 * n45 * n52 -
650 n15 * n24 * n42 * n53 + n14 * n25 * n42 * n53 + n15 * n22 * n44 * n53 -
651 n12 * n25 * n44 * n53 - n14 * n22 * n45 * n53 + n12 * n24 * n45 * n53 +
652 n15 * n23 * n42 * n54 - n13 * n25 * n42 * n54 - n15 * n22 * n43 * n54 +
653 n12 * n25 * n43 * n54 + n13 * n22 * n45 * n54 - n12 * n23 * n45 * n54 -
654 n14 * n23 * n42 * n55 + n13 * n24 * n42 * n55 + n14 * n22 * n43 * n55 -
655 n12 * n24 * n43 * n55 - n13 * n22 * n44 * n55 + n12 * n23 * n44 * n55 ) /
658 i[0][3] = ( -n15 * n24 * n33 * n52 + n14 * n25 * n33 * n52 + n15 * n23 * n34 * n52 -
659 n13 * n25 * n34 * n52 - n14 * n23 * n35 * n52 + n13 * n24 * n35 * n52 +
660 n15 * n24 * n32 * n53 - n14 * n25 * n32 * n53 - n15 * n22 * n34 * n53 +
661 n12 * n25 * n34 * n53 + n14 * n22 * n35 * n53 - n12 * n24 * n35 * n53 -
662 n15 * n23 * n32 * n54 + n13 * n25 * n32 * n54 + n15 * n22 * n33 * n54 -
663 n12 * n25 * n33 * n54 - n13 * n22 * n35 * n54 + n12 * n23 * n35 * n54 +
664 n14 * n23 * n32 * n55 - n13 * n24 * n32 * n55 - n14 * n22 * n33 * n55 +
665 n12 * n24 * n33 * n55 + n13 * n22 * n34 * n55 - n12 * n23 * n34 * n55 ) /
668 i[0][4] = ( n15 * n24 * n33 * n42 - n14 * n25 * n33 * n42 - n15 * n23 * n34 * n42 +
669 n13 * n25 * n34 * n42 + n14 * n23 * n35 * n42 - n13 * n24 * n35 * n42 -
670 n15 * n24 * n32 * n43 + n14 * n25 * n32 * n43 + n15 * n22 * n34 * n43 -
671 n12 * n25 * n34 * n43 - n14 * n22 * n35 * n43 + n12 * n24 * n35 * n43 +
672 n15 * n23 * n32 * n44 - n13 * n25 * n32 * n44 - n15 * n22 * n33 * n44 +
673 n12 * n25 * n33 * n44 + n13 * n22 * n35 * n44 - n12 * n23 * n35 * n44 -
674 n14 * n23 * n32 * n45 + n13 * n24 * n32 * n45 + n14 * n22 * n33 * n45 -
675 n12 * n24 * n33 * n45 - n13 * n22 * n34 * n45 + n12 * n23 * n34 * n45 ) /
678 double s0_prod = -0.07;
680 EvtComplex value0( 0., 0. );
681 EvtComplex value1( 0., 0. );
684 for (
int k = 0; k < 5; k++ )
686 double u1j_re =
real( i[0][k] );
687 double u1j_im =
imag( i[0][k] );
688 if ( u1j_re == 0. || u1j_im == 0. ) reject =
true;
691 for (
int pole_index = 0; pole_index < 5; pole_index++ )
693 EvtComplex
A = beta[pole_index] * g[pole_index][k];
694 value0 += ( i[0][k] *
A ) / ( ma[pole_index] * ma[pole_index] -
s );
698 value1 += i[0][k] * fprod[k];
702 value1 *= ( 1. - s0_prod ) / (
s - s0_prod );
704 if ( reject ==
true )
return EvtComplex( 9999., 9999. );
705 else return ( value0 + value1 );
712 double gammaR = 0.27;
714 if ( reso ==
"k0spim" ) mab2 = pow( ( p_k0l + p_pim ).
mass(), 2 );
715 else if ( reso ==
"k0spip" ) mab2 = pow( ( p_k0l + p_pip ).
mass(), 2 );
718 const double mD0 = 1.86483;
719 const double mKl = 0.49761;
720 const double mPi = 0.13957;
726 double _phiR = -1.9146;
727 double _phiF = 0.0017;
731 double mAB = sqrt( mab2 );
738 sqrt( ( ( ( mAB * mAB - mA * mA - mB * mB ) * ( mAB * mAB - mA * mA - mB * mB ) / 4.0 ) -
739 mA * mA * mB * mB ) /
744 sqrt( ( ( ( mR * mR - mA * mA - mB * mB ) * ( mR * mR - mA * mA - mB * mB ) / 4.0 ) -
745 mA * mA * mB * mB ) /
750 double g = gammaR * pow(
q / q0, power ) * ( mR / mAB ) * fR * fR;
751 EvtComplex propagator_relativistic_BreitWigner =
752 1. / ( mR * mR - mAB * mAB - EvtComplex( 0., mR * g ) );
755 double cot_deltaF = 1.0 / ( _a *
q ) + 0.5 * _r *
q;
756 double qcot_deltaF = 1.0 / _a + 0.5 * _r *
q *
q;
759 EvtComplex expi2deltaF = EvtComplex( qcot_deltaF,
q ) / EvtComplex( qcot_deltaF, -
q );
760 EvtComplex resonant_term_T =
761 _R * EvtComplex(
cos( _phiR + 2 * _phiF ),
sin( _phiR + 2 * _phiF ) ) *
762 propagator_relativistic_BreitWigner * mR * gammaR * mR / q0 * expi2deltaF;
765 EvtComplex non_resonant_term_F = _F * EvtComplex(
cos( _phiF ),
sin( _phiF ) ) *
766 (
cos( _phiF ) + cot_deltaF *
sin( _phiF ) ) * sqrt(
s ) /
767 EvtComplex( qcot_deltaF, -
q );
770 EvtComplex LASS_contribution = non_resonant_term_F + resonant_term_T;
772 return EvtComplex( A_r *
cos( Phi_r ), A_r *
sin( Phi_r ) ) * LASS_contribution;
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
Evt3Rank3C conj(const Evt3Rank3C &t2)
double imag(const EvtComplex &c)
double abs2(const EvtComplex &c)
ostream & report(Severity severity, const char *facility)
double sin(const BesAngle a)
double cos(const BesAngle a)
****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
***************************************************************************************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
static const double radToDegrees
void decay(EvtParticle *p)
void getName(std::string &name)
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)
void setProb(double prob)
const EvtVector4R & getP4() const
EvtParticle * getDaug(int i)
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)