74 width[0] = 5.0800e-02;
75 width[1] = 1.4041e-01;
76 width[2] = 2.6500e-01;
77 width[3] = 2.7000e-01;
78 width[4] = 5.0800e-02;
108 int GG[4][4] = { { 1, 0, 0, 0 }, { 0, -1, 0, 0 }, { 0, 0, -1, 0 }, { 0, 0, 0, -1 } };
109 for (
int i = 0; i < 4; i++ )
111 for (
int j = 0; j < 4; j++ ) { G[i][j] = GG[i][j]; }
155 double P1[4], P2[4], P3[4];
174 int g0[5] = { 1, 1, 1, 1, 12 };
175 int spin[5] = { 1, 0, 0, 0, 0 };
177 calEva( P1, P2, P3, mass, width, rho, phi, g0, spin, modetype, nstates, value );
183void EvtDsToKSKSpi::Com_Multi(
double a1[2],
double a2[2],
double res[2] ) {
184 res[0] = a1[0] * a2[0] - a1[1] * a2[1];
185 res[1] = a1[1] * a2[0] + a1[0] * a2[1];
187void EvtDsToKSKSpi::Com_Divide(
double a1[2],
double a2[2],
double res[2] ) {
188 double tmp = a2[0] * a2[0] + a2[1] * a2[1];
189 res[0] = ( a1[0] * a2[0] + a1[1] * a2[1] ) / tmp;
190 res[1] = ( a1[1] * a2[0] - a1[0] * a2[1] ) / tmp;
193double EvtDsToKSKSpi::SCADot(
double a1[4],
double a2[4] ) {
194 double _cal = a1[0] * a2[0] - a1[1] * a2[1] - a1[2] * a2[2] - a1[3] * a2[3];
197double EvtDsToKSKSpi::barrier(
int l,
double sa,
double sb,
double sc,
double r,
199 double q = ( sa + sb - sc ) * ( sa + sb - sc ) / ( 4 * sa ) - sb;
200 if (
q < 0 )
q = 1e-16;
205 double q0 = ( sa0 + sb - sc ) * ( sa0 + sb - sc ) / ( 4 * sa0 ) - sb;
206 if ( q0 < 0 ) q0 = 1e-16;
207 double z0 = q0 * r * r;
210 if ( l == 1 ) F = sqrt( ( 1 + z0 ) / ( 1 + z ) );
211 if ( l == 2 ) F = sqrt( ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ) );
215void EvtDsToKSKSpi::calt1(
double daug1[4],
double daug2[4],
double t1[4] ) {
218 for (
int i = 0; i < 4; i++ )
220 pa[i] = daug1[i] + daug2[i];
221 qa[i] = daug1[i] - daug2[i];
223 p = SCADot( pa, pa );
224 pq = SCADot( pa, qa );
226 for (
int i = 0; i < 4; i++ ) { t1[i] = qa[i] - tmp * pa[i]; }
228void EvtDsToKSKSpi::calt2(
double daug1[4],
double daug2[4],
double t2[4][4] ) {
231 calt1( daug1, daug2, t1 );
232 r = SCADot( t1, t1 ) / 3.0;
233 for (
int i = 0; i < 4; i++ ) { pa[i] = daug1[i] + daug2[i]; }
234 p = SCADot( pa, pa );
235 for (
int i = 0; i < 4; i++ )
237 for (
int j = 0; j < 4; j++ )
238 { t2[i][j] = t1[i] * t1[j] - r * ( G[i][j] - pa[i] * pa[j] / p ); }
243double EvtDsToKSKSpi::wid(
double mass2,
double mass,
double sa,
double sb,
double sc,
246 double m = sqrt( sa );
247 double tmp = sb - sc;
248 double tmp1 = sa + tmp;
249 double q = 0.25 * tmp1 * tmp1 / sa - sb;
250 if (
q < 0 )
q = 1e-16;
251 double tmp2 = mass2 + tmp;
252 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
253 if ( q0 < 0 ) q0 = 1e-16;
257 if ( l == 0 ) { widm = sqrt(
t ) * mass / m; }
258 else if ( l == 1 ) { widm =
t * sqrt(
t ) * mass / m * ( 1 + z0 ) / ( 1 + z ); }
260 { widm =
t *
t * sqrt(
t ) * mass / m * ( 9 + 3 * z0 + z0 * z0 ) / ( 9 + 3 * z + z * z ); }
263double EvtDsToKSKSpi::widl1(
double mass2,
double mass,
double sa,
double sb,
double sc,
266 double m = sqrt( sa );
267 double tmp = sb - sc;
268 double tmp1 = sa + tmp;
269 double q = 0.25 * tmp1 * tmp1 / sa - sb;
270 if (
q < 0 )
q = 1e-16;
271 double tmp2 = mass2 + tmp;
272 double q0 = 0.25 * tmp2 * tmp2 / mass2 - sb;
273 if ( q0 < 0 ) q0 = 1e-16;
276 double F = ( 1 + z0 ) / ( 1 + z );
278 widm =
t * sqrt(
t ) * mass / m * F;
281void EvtDsToKSKSpi::propagatorRBW(
double mass2,
double mass,
double width,
double sa,
282 double sb,
double sc,
double r2,
int l,
double prop[2] ) {
287 b[1] = -mass * width * wid( mass2, mass, sa, sb, sc, r2, l );
288 Com_Divide( a, b, prop );
291void EvtDsToKSKSpi::propagatorGS(
double mass2,
double mass,
double width,
double sa,
292 double sb,
double sc,
double r2,
double prop[2] ) {
294 double tmp = sb - sc;
295 double tmp1 = sa + tmp;
296 double q2 = 0.25 * tmp1 * tmp1 / sa - sb;
297 if ( q2 < 0 ) q2 = 1e-16;
299 double tmp2 = mass2 + tmp;
300 double q02 = 0.25 * tmp2 * tmp2 / mass2 - sb;
301 if ( q02 < 0 ) q02 = 1e-16;
303 double q = sqrt( q2 );
304 double q0 = sqrt( q02 );
305 double m = sqrt( sa );
306 double q03 = q0 * q02;
307 double tmp3 = log( mass + 2 * q0 ) + 1.2926305904;
309 double h = GS1 *
q / m * ( log( m + 2 *
q ) + 1.2926305904 );
310 double h0 = GS1 * q0 / mass * tmp3;
311 double dh = h0 * ( 0.125 / q02 - 0.5 / mass2 ) + GS3 / mass2;
312 double d = GS2 / q02 * tmp3 + GS3 * mass / q0 - GS4 * mass / q03;
313 double f = mass2 / q03 * ( q2 * ( h - h0 ) + ( mass2 - sa ) * q02 * dh );
315 a[0] = 1.0 + d * width / mass;
317 b[0] = mass2 - sa + width *
f;
318 b[1] = -mass * width * widl1( mass2, mass, sa, sb, sc, r2 );
319 Com_Divide( a, b, prop );
322void EvtDsToKSKSpi::KPiSLASS(
double sa,
double sb,
double sc,
double prop[2] ) {
323 const double m1430 = 1.441;
324 const double sa0 = 1.441 * 1.441;
325 const double w1430 = 0.193;
326 const double Lass1 = 0.25 / sa0;
327 double tmp = sb - sc;
328 double tmp1 = sa0 + tmp;
329 double q0 = Lass1 * tmp1 * tmp1 - sb;
330 if ( q0 < 0 ) q0 = 1e-16;
331 double tmp2 = sa + tmp;
332 double qs = 0.25 * tmp2 * tmp2 / sa - sb;
333 double q = sqrt( qs );
334 double width = w1430 *
q * m1430 / sqrt( sa * q0 );
335 double temp_R = atan( m1430 * width / ( sa0 - sa ) );
336 if ( temp_R < 0 ) temp_R += math_pi;
337 double deltaR = -109.7 + temp_R;
339 atan( 0.226 *
q / ( 2.0 - 3.819 * qs ) );
340 if ( temp_F < 0 ) temp_F += math_pi;
341 double deltaF = 0.1 + temp_F;
342 double deltaS = deltaR + 2.0 * deltaF;
343 double t1 = 0.96 *
sin( deltaF );
344 double t2 =
sin( deltaR );
346 CF[0] =
cos( deltaF );
347 CF[1] =
sin( deltaF );
348 CS[0] =
cos( deltaS );
349 CS[1] =
sin( deltaS );
350 prop[0] = t1 * CF[0] + t2 *
CS[0];
351 prop[1] = t1 * CF[1] + t2 *
CS[1];
354double EvtDsToKSKSpi::DDalitz(
double P1[4],
double P2[4],
double P3[4],
int Ang,
357 double temp_PDF, v_re;
360 double B[2], s1, s2, s3, sR, sD;
361 for (
int i = 0; i < 4; i++ )
363 pR[i] = P1[i] + P2[i];
364 pD[i] = pR[i] + P3[i];
366 s1 = SCADot( P1, P1 );
367 s2 = SCADot( P2, P2 );
368 s3 = SCADot( P3, P3 );
369 sR = SCADot( pR, pR );
370 sD = SCADot( pD, pD );
389 B[0] = barrier( 1, sR, s1, s2, 3.0, mass );
390 B[1] = barrier( 1, sD, sR, s3, 5.0, 1.9683 );
395 for (
int i = 0; i < 4; i++ ) { temp_PDF += t1[i] * T1[i] * G[i][i]; }
399 B[0] = barrier( 2, sR, s1, s2, 3.0, mass );
400 B[1] = barrier( 2, sD, sR, s3, 5.0, 1.9683 );
401 double t2[4][4], T2[4][4];
405 for (
int i = 0; i < 4; i++ )
407 for (
int j = 0; j < 4; j++ ) { temp_PDF += t2[i][j] * T2[j][i] * G[i][i] * G[j][j]; }
410 v_re = temp_PDF *
B[0] *
B[1];
414void EvtDsToKSKSpi::calEva(
double* Ks01,
double* Ks02,
double* Pic,
double* mass1,
415 double* width1,
double* amp,
double* phase,
int* g0,
int* spin,
416 int* modetype,
int nstates,
double& Result ) {
419 double P12[4], P23[4], P13[4];
420 double cof[2], amp_PDF[2], PDF[2];
421 double scpi, sks02, sks01;
422 double s12, s13, s23;
423 for (
int i = 0; i < 4; i++ )
425 P12[i] = Ks01[i] + Ks02[i];
426 P13[i] = Ks01[i] + Pic[i];
427 P23[i] = Ks02[i] + Pic[i];
429 scpi = SCADot( Pic, Pic );
430 sks01 = SCADot( Ks01, Ks01 );
431 sks02 = SCADot( Ks02, Ks02 );
432 s12 = SCADot( P12, P12 );
433 s13 = SCADot( P13, P13 );
434 s23 = SCADot( P23, P23 );
435 double pro[2], temp_PDF, amp_tmp[2], temp_PDF1, temp_PDF2, pro1[2], pro2[2];
443 for (
int i = 0; i < nstates; i++ )
447 mass1sq = mass1[i] * mass1[i];
448 cof[0] = amp[i] *
cos( phase[i] );
449 cof[1] = amp[i] *
sin( phase[i] );
463 if ( modetype[i] == 12 )
465 temp_PDF = DDalitz( Ks01, Ks02, Pic, spin[i], mass1[i] );
467 propagatorRBW( mass1sq, mass1[i], width1[i], s12, sks01, sks02, rRes, spin[i], pro );
468 if ( g0[i] == 12 ) KPiSLASS( s12, sks01, sks02, pro );
474 amp_tmp[0] = temp_PDF * pro[0];
475 amp_tmp[1] = temp_PDF * pro[1];
477 if ( modetype[i] == 13 )
479 temp_PDF1 = DDalitz( Ks01, Pic, Ks02, spin[i], mass1[i] );
481 propagatorRBW( mass1sq, mass1[i], width1[i], s13, sks01, scpi, rRes, spin[i], pro1 );
482 if ( g0[i] == 12 ) KPiSLASS( s13, sks01, scpi, pro1 );
488 temp_PDF2 = DDalitz( Ks02, Pic, Ks01, spin[i], mass1[i] );
490 propagatorRBW( mass1sq, mass1[i], width1[i], s23, sks02, scpi, rRes, spin[i], pro2 );
491 if ( g0[i] == 12 ) KPiSLASS( s23, sks02, scpi, pro2 );
497 amp_tmp[0] = ( temp_PDF1 * pro1[0] + temp_PDF2 * pro2[0] );
498 amp_tmp[1] = ( temp_PDF1 * pro1[1] + temp_PDF2 * pro2[1] );
500 Com_Multi( amp_tmp, cof, amp_PDF );
501 PDF[0] += amp_PDF[0];
502 PDF[1] += amp_PDF[1];
504 double value = PDF[0] * PDF[0] + PDF[1] * PDF[1];
506 if ( value <= 0 ) value = 1e-20;
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
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
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)
void decay(EvtParticle *p)
void getName(std::string &name)
const EvtVector4R & getP4() const
EvtParticle * getDaug(int i)
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)