107 int GG[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
108 for (
int i = 0; i < 4; i++ )
110 for (
int j = 0; j < 4; j++ ) { G[i][j] = GG[i][j]; }
123 double P1[4], P2[4], P3[4];
138 int g0[4] = { 1, 2, 2, 6 };
139 int spin[4] = { 2, 1, 1, 0 };
141 double r0[4] = { 3, 3, 3, 3 };
142 double r1[4] = { 5, 5, 5, 5 };
144 calEva( P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value, r0, r1 );
150void EvtDTopipi0pi0::Com_Multi(
double a1[2],
double a2[2],
double res[2] ) {
151 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
152 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
154void EvtDTopipi0pi0::Com_Divide(
double a1[2],
double a2[2],
double res[2] ) {
155 double tmp = a2[0] * a2[0] + a2[1] * a2[1];
156 res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / tmp;
157 res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / tmp;
160double EvtDTopipi0pi0::CalRho4pi( double_t
s ) {
161 if (
s >= 1. ) {
return sqrt( (
s - 16. * mass_Pion * mass_Pion ) /
s ); }
164 double_t s0 = 1.2274 + 0.00370909 / (
s *
s ) - ( 0.111203 ) / (
s)-6.39017 *
s +
165 16.8358 *
s *
s - 21.8845 *
s *
s *
s + 11.3153 *
s *
s *
s *
s;
166 double_t gam = s0 * sqrt( 1.0 - ( 16.0 * mass_Pion * mass_Pion ) );
173double EvtDTopipi0pi0::SCADot(
double a1[4],
double a2[4] ) {
174 double _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
177double EvtDTopipi0pi0::barrier(
int l,
double sa,
double sb,
double sc,
double r,
179 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
180 if (
q < 0 )
q = 1e-16;
185 double q0 = ( sa0 + sb - sc ) * ( sa0 + sb - sc ) / ( 4 * sa0 ) - sb;
186 if ( q0 < 0 ) q0 = 1e-16;
187 double z0 = q0 * r * r;
190 if ( l == 1 ) F = sqrt( ( 1 + z0 ) / ( 1 + z ) );
191 if ( l == 2 ) F = sqrt( ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ) );
195void EvtDTopipi0pi0::calt1(
double daug1[4],
double daug2[4],
double t1[4] ) {
198 for (
int i = 0; i < 4; i++ )
200 pa[i] = daug1[i] + daug2[i];
201 qa[i] = daug1[i] - daug2[i];
203 p = SCADot( pa, pa );
204 pq = SCADot( pa, qa );
206 for (
int i = 0; i < 4; i++ ) { t1[i] = qa[i] - tmp * pa[i]; }
208void EvtDTopipi0pi0::calt2(
double daug1[4],
double daug2[4],
double t2[4][4] ) {
211 calt1( daug1, daug2, t1 );
212 r = SCADot( t1, t1 ) / 3.0;
213 for (
int i = 0; i < 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
214 p = SCADot( pa, pa );
215 for (
int i = 0; i < 4; i++ )
217 for (
int j = 0; j < 4; j++ )
218 { t2[i][j] = t1[i] * t1[j] - r * ( G[i][j] - pa[i] * pa[j] / p ); }
222void EvtDTopipi0pi0::propagator(
double mass2,
double mass,
double width,
double sx,
228 b[1] = -mass * width;
229 Com_Divide( a, b, prop );
231double EvtDTopipi0pi0::wid(
double mass2,
double mass,
double sa,
double sb,
double sc,
234 double m = sqrt( sa );
235 double tmp = sb - sc;
236 double tmp1 = sa + tmp;
237 double q = 0.25 * tmp1 * tmp1 / sa - sb;
238 if (
q < 0 )
q = 1e-16;
239 double tmp2 = mass2 + tmp;
240 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
241 if ( q0 < 0 ) q0 = 1e-16;
245 if ( l == 0 ) { widm = sqrt(
t ) * mass / m; }
246 else if ( l == 1 ) { widm =
t * sqrt(
t ) * mass / m * ( 1 + z0 ) / ( 1 + z ); }
248 { widm =
t *
t * sqrt(
t ) * mass / m * ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ); }
251double EvtDTopipi0pi0::widl1(
double mass2,
double mass,
double sa,
double sb,
double sc,
254 double m = sqrt( sa );
255 double tmp = sb - sc;
256 double tmp1 = sa + tmp;
257 double q = 0.25 * tmp1 * tmp1 / sa - sb;
258 if (
q < 0 )
q = 1e-16;
259 double tmp2 = mass2 + tmp;
260 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
261 if ( q0 < 0 ) q0 = 1e-16;
264 double F = ( 1 + z0 ) / ( 1 + z );
266 widm =
t * sqrt(
t ) * mass / m * F;
269void EvtDTopipi0pi0::propagatorRBW(
double mass2,
double mass,
double width,
double sa,
270 double sb,
double sc,
double r2,
int l,
double prop[2] ) {
276 b[1] = -mass * width * wid( mass2, mass, sa, sb, sc, r2, l );
277 Com_Divide( a, b, prop );
280void EvtDTopipi0pi0::propagatorFlatte(
double mass,
double width,
double sa,
double sb,
281 double sc,
int r,
double prop[2] ) {
283 double rhoab[2], rhoKsK[2];
284 q = 0.25 * ( sa + sb - sc ) * ( sa + sb - sc ) / sa - sb;
285 if ( r == 0 ) qKsK = 0.25 * sa - 0.49368 * 0.49368;
288 0.25 * ( sa + 3.899750596e-03 ) * ( sa + 3.899750596e-03 ) / sa - 0.497614 * 0.497614;
291 rhoab[0] = 2 * sqrt(
q / sa );
297 rhoab[1] = 2 * sqrt( -
q / sa );
301 rhoKsK[0] = 2 * sqrt( qKsK / sa );
307 rhoKsK[1] = 2 * sqrt( -qKsK / sa );
312 b[0] = mass * mass - sa + 0.341 * rhoab[1] + 0.892 * 0.341 * rhoKsK[1];
313 b[1] = -( 0.341 * rhoab[0] + 0.892 * 0.341 * rhoKsK[0] );
314 Com_Divide( a, b, prop );
317void EvtDTopipi0pi0::rhoab(
double sa,
double sb,
double sc,
double res[2] ) {
318 double tmp = sa + sb - sc;
319 double q = 0.25 * tmp * tmp / sa - sb;
322 res[0] = 2.0 * sqrt(
q / sa );
328 res[1] = 2.0 * sqrt( -
q / sa );
331void EvtDTopipi0pi0::rho4Pi(
double sa,
double res[2] ) {
332 double temp = 1.0 - 0.3116765584 / sa;
335 res[0] = sqrt( temp ) / ( 1.0 +
exp( 9.8 - 3.5 * sa ) );
341 res[1] = sqrt( -temp ) / ( 1.0 +
exp( 9.8 - 3.5 * sa ) );
345void EvtDTopipi0pi0::propagatorsigma500(
double sa,
double sb,
double sc,
double prop[2] ) {
346 double f = 0.5843 + 1.6663 * sa;
347 const double M = 0.9264;
348 const double mass2 = 0.85821696;
349 const double mpi2d2 = 0.00973989245;
350 double g1 =
f * ( sa - mpi2d2 ) / ( mass2 - mpi2d2 ) *
exp( ( mass2 - sa ) / 1.082 );
351 double rho1s[2], rho1M[2], rho2s[2], rho2M[2], rho1[2], rho2[2];
352 rhoab( sa, sb, sc, rho1s );
353 rhoab( mass2, sb, sc, rho1M );
355 rho4Pi( mass2, rho2M );
356 Com_Divide( rho1s, rho1M, rho1 );
357 Com_Divide( rho2s, rho2M, rho2 );
361 b[0] = mass2 - sa + M * (
g1 * rho1[1] + 0.0024 * rho2[1] );
362 b[1] = -M * (
g1 * rho1[0] + 0.0024 * rho2[0] );
363 Com_Divide( a, b, prop );
365void EvtDTopipi0pi0::Flatte_rhoab(
double sa,
double sb,
double rho[2] ) {
366 double q = 1.0 - ( 4 * sb / sa );
380void EvtDTopipi0pi0::propagator980(
double mass,
double sx,
double* sb,
double prop[2] ) {
382 double gpipi1 = 2.0 / 3.0 * 0.165;
383 double gpipi2 = 1.0 / 3.0 * 0.165;
384 double gKK1 = 0.5 * 0.69465;
385 double gKK2 = 0.5 * 0.69465;
387 double unit[2] = { 1.0 };
388 double ci[2] = { 0, 1 };
390 Flatte_rhoab( sx, sb[0], rho1 );
392 Flatte_rhoab( sx, sb[1], rho2 );
394 Flatte_rhoab( sx, sb[2], rho3 );
396 Flatte_rhoab( sx, sb[3], rho4 );
398 double tmp1[2] = { gpipi1, 0 };
400 double tmp2[2] = { gpipi2, 0 };
402 double tmp3[2] = { gKK1, 0 };
404 double tmp4[2] = { gKK2, 0 };
407 Com_Multi( tmp1, rho1, tmp11 );
408 Com_Multi( tmp2, rho2, tmp22 );
409 Com_Multi( tmp3, rho3, tmp33 );
410 Com_Multi( tmp4, rho4, tmp44 );
412 double tmp5[2] = { tmp11[0] + tmp22[0] + tmp33[0] + tmp44[0],
413 tmp11[1] + tmp22[1] + tmp33[1] + tmp44[1] };
415 Com_Multi( tmp5, ci, tmp51 );
416 double tmp6[2] = { mass * mass - sx - tmp51[0], -1.0 * tmp51[1] };
418 Com_Divide(
unit, tmp6, prop );
422void EvtDTopipi0pi0::propagatorGS(
double mass2,
double mass,
double width,
double sa,
423 double sb,
double sc,
double r2,
double prop[2] ) {
425 double tmp = sb - sc;
426 double tmp1 = sa + tmp;
427 double q2 = 0.25 * tmp1 * tmp1 / sa - sb;
428 if ( q2 < 0 ) q2 = 1e-16;
430 double tmp2 = mass2 + tmp;
431 double q02 = 0.25 * tmp2 * tmp2 / mass2 - sb;
432 if ( q02 < 0 ) q02 = 1e-16;
434 double q = sqrt( q2 );
435 double q0 = sqrt( q02 );
436 double m = sqrt( sa );
437 double q03 = q0 * q02;
438 double tmp3 = log( mass + 2 * q0 ) + 1.2926305904;
440 double h = GS1 *
q / m * ( log( m + 2 *
q ) + 1.2926305904 );
441 double h0 = GS1 * q0 / mass * tmp3;
442 double dh = h0 * ( 0.125 / q02 - 0.5 / mass2 ) + GS3 / mass2;
443 double d = GS2 / q02 * tmp3 + GS3 * mass / q0 - GS4 * mass / q03;
444 double f = mass2 / q03 * ( q2 * ( h - h0 ) + ( mass2 - sa ) * q02 * dh );
446 a[0] = 1.0 + d * width / mass;
448 b[0] = mass2 - sa + width *
f;
449 b[1] = -mass * width * widl1( mass2, mass, sa, sb, sc, r2 );
450 Com_Divide( a, b, prop );
453double EvtDTopipi0pi0::DDalitz(
double P1[4],
double P2[4],
double P3[4],
int Ang,
double mass,
454 double rRES0,
double rRES1 ) {
456 double temp_PDF, v_re;
459 double B[2], s1, s2, s3, sR, sD;
460 for (
int i = 0; i < 4; i++ )
462 pR[i] = P1[i] + P2[i];
463 pD[i] = pR[i] + P3[i];
465 s1 = SCADot( P1, P1 );
466 s2 = SCADot( P2, P2 );
467 s3 = SCADot( P3, P3 );
468 sR = SCADot( pR, pR );
469 sD = SCADot( pD, pD );
471 for (
int i = 0; i != 4; i++ )
473 for (
int j = 0; j != 4; j++ )
477 if ( i == 0 ) G[i][j] = 1;
491 B[0] = barrier( 1, sR, s1, s2, rRES0, mass );
492 B[1] = barrier( 1, sD, sR, s3, rRES1, 1.86966 );
497 for (
int i = 0; i < 4; i++ ) { temp_PDF += t1[i] * T1[i] * G[i][i]; }
501 B[0] = barrier( 2, sR, s1, s2, rRES0, mass );
502 B[1] = barrier( 2, sD, sR, s3, rRES1, 1.86966 );
503 double t2[4][4], T2[4][4];
507 for (
int i = 0; i < 4; i++ )
509 for (
int j = 0; j < 4; j++ ) { temp_PDF += t2[i][j] * T2[j][i] * G[i][i] * G[j][j]; }
512 v_re = temp_PDF *
B[0] *
B[1];
516TComplex EvtDTopipi0pi0::ResonanceSkm( Double_t&
m2 ) {
517 Double_t g11 = 0.22889, g12 = -0.55377, g13 = 0.00, g14 = -0.39899, g15 = -0.34639;
518 Double_t g21 = 0.94128, g22 = 0.55095, g23 = 0.00, g24 = 0.39065, g25 = 0.31503;
519 Double_t g31 = 0.36856, g32 = 0.23888, g33 = 0.55639, g34 = 0.18340, g35 = 0.18681;
520 Double_t g41 = 0.33650, g42 = 0.40907, g43 = 0.85679, g44 = 0.19906, g45 = -0.00984;
521 Double_t g51 = 0.18171, g52 = -0.17558, g53 = -0.79658, g54 = -0.00355, g55 = 0.22358;
523 Double_t fs11 = 0.23399, fs12 = 0.15044, fs13 = -0.20545, fs14 = 0.32825, fs15 = 0.35412;
524 Double_t m_1 = 0.65100, m_2 = 1.20360, m_3 = 1.55817, m_4 = 1.21000, m_5 = 1.82206;
525 Double_t ss0 = -3.92637;
526 Double_t sp0 = -0.07;
527 double_t sA0 = -0.15;
530 double_t km11 = ( g11 * g11 / ( m_1 * m_1 -
m2 ) + g21 * g21 / ( m_2 * m_2 -
m2 ) +
531 g31 * g31 / ( m_3 * m_3 -
m2 ) + g41 * g41 / ( m_4 * m_4 -
m2 ) +
532 g51 * g51 / ( m_5 * m_5 -
m2 ) + fs11 * ( 1 - ss0 ) / (
m2 - ss0 ) ) *
533 ( 1 - sA0 ) / (
m2 - sA0 ) * (
m2 - sA * mass_Pion * mass_Pion * 0.5 );
534 double_t km12 = ( g11 * g12 / ( m_1 * m_1 -
m2 ) + g21 * g22 / ( m_2 * m_2 -
m2 ) +
535 g31 * g32 / ( m_3 * m_3 -
m2 ) + g41 * g42 / ( m_4 * m_4 -
m2 ) +
536 g51 * g52 / ( m_5 * m_5 -
m2 ) + fs12 * ( 1 - ss0 ) / (
m2 - ss0 ) ) *
537 ( 1 - sA0 ) / (
m2 - sA0 ) * (
m2 - sA * mass_Pion * mass_Pion * 0.5 );
538 double_t km13 = ( g11 * g13 / ( m_1 * m_1 -
m2 ) + g21 * g23 / ( m_2 * m_2 -
m2 ) +
539 g31 * g33 / ( m_3 * m_3 -
m2 ) + g41 * g43 / ( m_4 * m_4 -
m2 ) +
540 g51 * g53 / ( m_5 * m_5 -
m2 ) + fs13 * ( 1 - ss0 ) / (
m2 - ss0 ) ) *
541 ( 1 - sA0 ) / (
m2 - sA0 ) * (
m2 - sA * mass_Pion * mass_Pion * 0.5 );
542 double_t km14 = ( g11 * g14 / ( m_1 * m_1 -
m2 ) + g21 * g24 / ( m_2 * m_2 -
m2 ) +
543 g31 * g34 / ( m_3 * m_3 -
m2 ) + g41 * g44 / ( m_4 * m_4 -
m2 ) +
544 g51 * g54 / ( m_5 * m_5 -
m2 ) + fs14 * ( 1 - ss0 ) / (
m2 - ss0 ) ) *
545 ( 1 - sA0 ) / (
m2 - sA0 ) * (
m2 - sA * mass_Pion * mass_Pion * 0.5 );
546 double_t km15 = ( g11 * g15 / ( m_1 * m_1 -
m2 ) + g21 * g25 / ( m_2 * m_2 -
m2 ) +
547 g31 * g35 / ( m_3 * m_3 -
m2 ) + g41 * g45 / ( m_4 * m_4 -
m2 ) +
548 g51 * g55 / ( m_5 * m_5 -
m2 ) + fs15 * ( 1 - ss0 ) / (
m2 - ss0 ) ) *
549 ( 1 - sA0 ) / (
m2 - sA0 ) * (
m2 - sA * mass_Pion * mass_Pion * 0.5 );
551 double_t km21 = km12, km31 = km13, km41 = km14, km51 = km15;
552 double_t km23 = ( g12 * g13 / ( m_1 * m_1 -
m2 ) + g22 * g23 / ( m_2 * m_2 -
m2 ) +
553 g32 * g33 / ( m_3 * m_3 -
m2 ) + g42 * g43 / ( m_4 * m_4 -
m2 ) +
554 g52 * g53 / ( m_5 * m_5 -
m2 ) ) *
555 ( 1 - sA0 ) / (
m2 - sA0 ) * (
m2 - sA * mass_Pion * mass_Pion * 0.5 );
556 double_t km24 = ( g12 * g14 / ( m_1 * m_1 -
m2 ) + g22 * g24 / ( m_2 * m_2 -
m2 ) +
557 g32 * g34 / ( m_3 * m_3 -
m2 ) + g42 * g44 / ( m_4 * m_4 -
m2 ) +
558 g52 * g54 / ( m_5 * m_5 -
m2 ) ) *
559 ( 1 - sA0 ) / (
m2 - sA0 ) * (
m2 - sA * mass_Pion * mass_Pion * 0.5 );
560 double_t km25 = ( g12 * g15 / ( m_1 * m_1 -
m2 ) + g22 * g25 / ( m_2 * m_2 -
m2 ) +
561 g32 * g35 / ( m_3 * m_3 -
m2 ) + g42 * g45 / ( m_4 * m_4 -
m2 ) +
562 g52 * g55 / ( m_5 * m_5 -
m2 ) ) *
563 ( 1 - sA0 ) / (
m2 - sA0 ) * (
m2 - sA * mass_Pion * mass_Pion * 0.5 );
564 double_t km32 = km23, km42 = km24, km52 = km25;
565 double_t km34 = ( g13 * g14 / ( m_1 * m_1 -
m2 ) + g23 * g24 / ( m_2 * m_2 -
m2 ) +
566 g33 * g34 / ( m_3 * m_3 -
m2 ) + g43 * g44 / ( m_4 * m_4 -
m2 ) +
567 g53 * g54 / ( m_5 * m_5 -
m2 ) ) *
568 ( 1 - sA0 ) / (
m2 - sA0 ) * (
m2 - sA * mass_Pion * mass_Pion * 0.5 );
569 double_t km35 = ( g13 * g15 / ( m_1 * m_1 -
m2 ) + g23 * g25 / ( m_2 * m_2 -
m2 ) +
570 g33 * g35 / ( m_3 * m_3 -
m2 ) + g43 * g45 / ( m_4 * m_4 -
m2 ) +
571 g53 * g55 / ( m_5 * m_5 -
m2 ) ) *
572 ( 1 - sA0 ) / (
m2 - sA0 ) * (
m2 - sA * mass_Pion * mass_Pion * 0.5 );
573 double_t km43 = km34, km53 = km35;
574 double_t km45 = ( g14 * g15 / ( m_1 * m_1 -
m2 ) + g24 * g25 / ( m_2 * m_2 -
m2 ) +
575 g34 * g35 / ( m_3 * m_3 -
m2 ) + g44 * g45 / ( m_4 * m_4 -
m2 ) +
576 g54 * g55 / ( m_5 * m_5 -
m2 ) ) *
577 ( 1 - sA0 ) / (
m2 - sA0 ) * (
m2 - sA * mass_Pion * mass_Pion * 0.5 );
578 double_t km54 = km45;
579 double_t km22 = ( g12 * g12 / ( m_1 * m_1 -
m2 ) + g22 * g22 / ( m_2 * m_2 -
m2 ) +
580 g32 * g32 / ( m_3 * m_3 -
m2 ) + g42 * g42 / ( m_4 * m_4 -
m2 ) +
581 g52 * g52 / ( m_5 * m_5 -
m2 ) ) *
582 ( 1 - sA0 ) / (
m2 - sA0 ) * (
m2 - sA * mass_Pion * mass_Pion * 0.5 );
583 double_t km33 = ( g13 * g13 / ( m_1 * m_1 -
m2 ) + g23 * g23 / ( m_2 * m_2 -
m2 ) +
584 g33 * g33 / ( m_3 * m_3 -
m2 ) + g43 * g43 / ( m_4 * m_4 -
m2 ) +
585 g53 * g53 / ( m_5 * m_5 -
m2 ) ) *
586 ( 1 - sA0 ) / (
m2 - sA0 ) * (
m2 - sA * mass_Pion * mass_Pion * 0.5 );
587 double_t km44 = ( g14 * g14 / ( m_1 * m_1 -
m2 ) + g24 * g24 / ( m_2 * m_2 -
m2 ) +
588 g34 * g34 / ( m_3 * m_3 -
m2 ) + g44 * g44 / ( m_4 * m_4 -
m2 ) +
589 g54 * g54 / ( m_5 * m_5 -
m2 ) ) *
590 ( 1 - sA0 ) / (
m2 - sA0 ) * (
m2 - sA * mass_Pion * mass_Pion * 0.5 );
591 double_t km55 = ( g15 * g15 / ( m_1 * m_1 -
m2 ) + g25 * g25 / ( m_2 * m_2 -
m2 ) +
592 g35 * g35 / ( m_3 * m_3 -
m2 ) + g45 * g45 / ( m_4 * m_4 -
m2 ) +
593 g55 * g55 / ( m_5 * m_5 -
m2 ) ) *
594 ( 1 - sA0 ) / (
m2 - sA0 ) * (
m2 - sA * mass_Pion * mass_Pion * 0.5 );
596 TComplex fp11( 1.41620 *
cos( 5.90540 ), 1.41620 *
sin( 5.90540 ) );
597 TComplex fp12( 15.3640 *
cos( 2.22510 ), 15.3640 *
sin( 2.22510 ) );
598 TComplex fp13( 8.2297 *
cos( 2.30550 ), 8.2297 *
sin( 2.30550 ) );
602 TComplex beta1( 3.70140 *
cos( 5.0818 ), 3.70140 *
sin( 5.0818 ) );
603 TComplex beta2( 4.04470 *
cos( 0.13605 ), 4.04470 *
sin( 0.13605 ) );
604 TComplex beta3( 3.4591 *
cos( -5.6845 ), 3.4591 *
sin( -5.6845 ) );
608 TComplex P1 = fp11 * ( 1 - sp0 ) / (
m2 - sp0 ) + beta1 * g11 / ( m_1 * m_1 -
m2 ) +
609 beta2 * g21 / ( m_2 * m_2 -
m2 ) + beta3 * g31 / ( m_3 * m_3 -
m2 ) +
610 beta4 * g41 / ( m_4 * m_4 -
m2 ) + beta5 * g51 / ( m_5 * m_5 -
m2 );
611 TComplex P2 = fp12 * ( 1 - sp0 ) / (
m2 - sp0 ) + beta1 * g12 / ( m_1 * m_1 -
m2 ) +
612 beta2 * g22 / ( m_2 * m_2 -
m2 ) + beta3 * g32 / ( m_3 * m_3 -
m2 ) +
613 beta4 * g42 / ( m_4 * m_4 -
m2 ) + beta5 * g52 / ( m_5 * m_5 -
m2 );
614 TComplex P3 = fp13 * ( 1 - sp0 ) / (
m2 - sp0 ) + beta1 * g13 / ( m_1 * m_1 -
m2 ) +
615 beta2 * g23 / ( m_2 * m_2 -
m2 ) + beta3 * g33 / ( m_3 * m_3 -
m2 ) +
616 beta4 * g43 / ( m_4 * m_4 -
m2 ) + beta5 * g53 / ( m_5 * m_5 -
m2 );
617 TComplex P4 = fp14 * ( 1 - sp0 ) / (
m2 - sp0 ) + beta1 * g14 / ( m_1 * m_1 -
m2 ) +
618 beta2 * g24 / ( m_2 * m_2 -
m2 ) + beta3 * g34 / ( m_3 * m_3 -
m2 ) +
619 beta4 * g44 / ( m_4 * m_4 -
m2 ) + beta5 * g54 / ( m_5 * m_5 -
m2 );
620 TComplex P5 = fp15 * ( 1 - sp0 ) / (
m2 - sp0 ) + beta1 * g15 / ( m_1 * m_1 -
m2 ) +
621 beta2 * g25 / ( m_2 * m_2 -
m2 ) + beta3 * g35 / ( m_3 * m_3 -
m2 ) +
622 beta4 * g45 / ( m_4 * m_4 -
m2 ) + beta5 * g55 / ( m_5 * m_5 -
m2 );
624 TMatrixD MI( 5, 5 ), MA( 5, 5 ), MA_invt( 5, 5 ), MB( 5, 5 ), KM( 5, 5 ), RhoA( 5, 5 ),
625 RhoB( 5, 5 ), MRe( 5, 5 ), MIm( 5, 5 );
629 TMatrixDRow( KM, 0 )( 0 ) = km11;
630 TMatrixDRow( KM, 1 )( 0 ) = km21;
631 TMatrixDRow( KM, 2 )( 0 ) = km31;
632 TMatrixDRow( KM, 3 )( 0 ) = km41;
633 TMatrixDRow( KM, 4 )( 0 ) = km51;
634 TMatrixDRow( KM, 0 )( 1 ) = km12;
635 TMatrixDRow( KM, 1 )( 1 ) = km22;
636 TMatrixDRow( KM, 2 )( 1 ) = km32;
637 TMatrixDRow( KM, 3 )( 1 ) = km42;
638 TMatrixDRow( KM, 4 )( 1 ) = km52;
639 TMatrixDRow( KM, 0 )( 2 ) = km13;
640 TMatrixDRow( KM, 1 )( 2 ) = km23;
641 TMatrixDRow( KM, 2 )( 2 ) = km33;
642 TMatrixDRow( KM, 3 )( 2 ) = km43;
643 TMatrixDRow( KM, 4 )( 2 ) = km53;
644 TMatrixDRow( KM, 0 )( 3 ) = km14;
645 TMatrixDRow( KM, 1 )( 3 ) = km24;
646 TMatrixDRow( KM, 2 )( 3 ) = km34;
647 TMatrixDRow( KM, 3 )( 3 ) = km44;
648 TMatrixDRow( KM, 4 )( 3 ) = km54;
649 TMatrixDRow( KM, 0 )( 4 ) = km15;
650 TMatrixDRow( KM, 1 )( 4 ) = km25;
651 TMatrixDRow( KM, 2 )( 4 ) = km35;
652 TMatrixDRow( KM, 3 )( 4 ) = km45;
653 TMatrixDRow( KM, 4 )( 4 ) = km55;
655 if ( ( 1.0 - 4.0 * mass_Pion * mass_Pion /
m2 ) > 0 )
657 TMatrixDRow( RhoA, 0 )( 0 ) = sqrt( 1.0 - 4.0 * mass_Pion * mass_Pion /
m2 );
658 TMatrixDRow( RhoB, 0 )( 0 ) = 0.0;
662 TMatrixDRow( RhoA, 0 )( 0 ) = 0.0;
663 TMatrixDRow( RhoB, 0 )( 0 ) = sqrt( 4.0 * mass_Pion * mass_Pion /
m2 - 1.0 );
666 if ( ( 1.0 - 4.0 * mass_Kaon * mass_Kaon /
m2 ) > 0 )
668 TMatrixDRow( RhoA, 1 )( 1 ) = sqrt( 1.0 - 4.0 * mass_Kaon * mass_Kaon /
m2 );
669 TMatrixDRow( RhoB, 1 )( 1 ) = 0.0;
673 TMatrixDRow( RhoA, 1 )( 1 ) = 0.0;
674 TMatrixDRow( RhoB, 1 )( 1 ) = sqrt( 4.0 * mass_Kaon * mass_Kaon /
m2 - 1.0 );
677 TMatrixDRow( RhoA, 2 )( 2 ) = CalRho4pi(
m2 );
678 TMatrixDRow( RhoB, 2 )( 2 ) = 0.0;
679 if ( ( 1.0 - 4.0 * mass_Eta * mass_Eta /
m2 ) > 0 )
681 TMatrixDRow( RhoA, 3 )( 3 ) = sqrt( 1.0 - 4.0 * mass_Eta * mass_Eta /
m2 );
682 TMatrixDRow( RhoB, 3 )( 3 ) = 0.0;
686 TMatrixDRow( RhoA, 3 )( 3 ) = 0.0;
687 TMatrixDRow( RhoB, 3 )( 3 ) = sqrt( 4.0 * mass_Eta * mass_Eta /
m2 - 1.0 );
690 if ( ( 1.0 - ( mass_Eta + mass_Etap ) * ( mass_Eta + mass_Etap ) /
m2 ) > 0 )
692 TMatrixDRow( RhoA, 4 )( 4 ) =
693 sqrt( 1.0 - ( mass_Eta + mass_Etap ) * ( mass_Eta + mass_Etap ) /
m2 );
694 TMatrixDRow( RhoB, 4 )( 4 ) = 0.0;
698 TMatrixDRow( RhoA, 4 )( 4 ) = 0.0;
699 TMatrixDRow( RhoB, 4 )( 4 ) =
700 sqrt( ( mass_Eta + mass_Etap ) * ( mass_Eta + mass_Etap ) /
m2 - 1.0 );
704 MB = -1.0 * KM * RhoA;
709 MRe = MA + MB * MA_invt * MB;
711 MIm = MA_invt * MB * MRe;
717 amp = ( ciR * TMatrixDRow( MRe, 0 )( 0 ) - ciM * TMatrixDRow( MIm, 0 )( 0 ) ) * P1 +
718 ( ciR * TMatrixDRow( MRe, 0 )( 1 ) - ciM * TMatrixDRow( MIm, 0 )( 1 ) ) * P2 +
719 ( ciR * TMatrixDRow( MRe, 0 )( 2 ) - ciM * TMatrixDRow( MIm, 0 )( 2 ) ) * P3 +
720 ( ciR * TMatrixDRow( MRe, 0 )( 3 ) - ciM * TMatrixDRow( MIm, 0 )( 3 ) ) * P4 +
721 ( ciR * TMatrixDRow( MRe, 0 )( 4 ) - ciM * TMatrixDRow( MIm, 0 )( 4 ) ) * P5;
726void EvtDTopipi0pi0::calEva(
double* Pic,
double* Pi01,
double* Pi02,
double* mass1,
727 double* width1,
double* amp,
double* phase,
int* g0,
int* spin,
728 int* modetype,
int nstates,
double& Result,
double* r0,
730 double P12[4], P23[4], P13[4];
732 double cof[2], amp_PDF[2], PDF[2];
733 double spi01, spi02, scpi;
735 double s12, s13, s23;
737 double Realpipis, Imagpipis;
738 for (
int i = 0; i < 4; i++ )
740 P23[i] = Pi01[i] + Pi02[i];
741 P12[i] = Pic[i] + Pi01[i];
742 P13[i] = Pic[i] + Pi02[i];
744 scpi = SCADot( Pic, Pic );
745 spi01 = SCADot( Pi01, Pi01 );
746 spi02 = SCADot( Pi02, Pi02 );
747 s23 = SCADot( P23, P23 );
748 s12 = SCADot( P12, P12 );
749 s13 = SCADot( P13, P13 );
751 double spi012[4] = { mass_Pion * mass_Pion, spi02, mass_Kaon * mass_Kaon, mk0 * mk0 };
753 double pro[2], temp_PDF, pro1[2], temp_PDF1, pro2[2], temp_PDF2, amp_tmp[2];
761 for (
int i = 0; i < nstates; i++ )
765 mass1sq = mass1[i] * mass1[i];
766 cof[0] = amp[i] *
cos( phase[i] );
767 cof[1] = amp[i] *
sin( phase[i] );
769 if ( modetype[i] == 23 )
771 temp_PDF = DDalitz( Pi01, Pi02, Pic, spin[i], mass1[i], r0[i], r1[i] );
772 if ( g0[i] == 2 ) propagatorsigma500( s23, spi01, spi02, pro );
774 propagatorRBW( mass1sq, mass1[i], width1[i], s23, spi01, spi02, r0[i] * r0[i], spin[i],
776 if ( g0[i] == 3 ) propagator980( mass1[i], s23, spi012, pro );
779 tmpswave = ResonanceSkm( s23 );
780 Realpipis = tmpswave.Re();
781 Imagpipis = tmpswave.Im();
782 amp_tmp[0] = temp_PDF * Realpipis;
783 amp_tmp[1] = temp_PDF * Imagpipis;
790 if ( g0[i] != 6 ) amp_tmp[0] = temp_PDF * pro[0];
791 if ( g0[i] != 6 ) amp_tmp[1] = temp_PDF * pro[1];
793 if ( modetype[i] == 12 )
795 temp_PDF1 = DDalitz( Pic, Pi01, Pi02, spin[i], mass1[i], r0[i], r1[i] );
797 propagatorRBW( mass1sq, mass1[i], width1[i], s12, scpi, spi01, r0[i] * r0[i], spin[i],
800 propagatorGS( mass1sq, mass1[i], width1[i], s12, scpi, spi01, r0[i] * r0[i], pro1 );
807 temp_PDF2 = DDalitz( Pic, Pi02, Pi01, spin[i], mass1[i], r0[i], r1[i] );
809 propagatorRBW( mass1sq, mass1[i], width1[i], s13, scpi, spi02, r0[i] * r0[i], spin[i],
812 propagatorGS( mass1sq, mass1[i], width1[i], s13, scpi, spi02, r0[i] * r0[i], pro2 );
818 amp_tmp[0] = temp_PDF1 * pro1[0] + temp_PDF2 * pro2[0];
820 amp_tmp[1] = temp_PDF1 * pro1[1] + temp_PDF2 * pro2[1];
822 Com_Multi( amp_tmp, cof, amp_PDF );
823 PDF[0] += amp_PDF[0];
824 PDF[1] += amp_PDF[1];
826 double value = PDF[0] * PDF[0] + PDF[1] * PDF[1];
827 if ( value <= 0 ) value = 1e-20;
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
complex< double > TComplex
EvtComplex exp(const EvtComplex &c)
*******INTEGER m_nBinMax INTEGER m_NdiMax !No of bins in histogram for cell exploration division $ !Last vertex $ !Last active cell $ !Last cell in buffer $ !No of sampling when dividing cell $ !No of function total $ !Flag for random ceel for $ !Flag for type of for WtMax $ !Flag which decides whether vertices are included in the sampling $ entire domain is hyp !Maximum effective eevents per saves r n generator level $ !Flag for chat level in !Latex Output unit
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
virtual ~EvtDTopipi0pi0()
void decay(EvtParticle *p)
void getName(std::string &name)
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)
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)