51 width0_1430_S = 0.270;
85 Pi = atan2( 0.0, -1.0 );
87 root2d3 = sqrt( 2. / 3 );
88 root1d2 = sqrt( 0.5 );
89 root3d2 = sqrt( 1.5 );
126 if ( pid == -321 ) 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 );
137 double& cosV,
double& cosL,
double& chi ) {
145 boost.
set( vp4_KPi.
get( 0 ), -vp4_KPi.
get( 1 ), -vp4_KPi.
get( 2 ), -vp4_KPi.
get( 3 ) );
147 cosV = vp4_Kp.
dot( vp4_KPi ) / ( vp4_Kp.
d3mag() * vp4_KPi.
d3mag() );
149 boost.
set( vp4_W.
get( 0 ), -vp4_W.
get( 1 ), -vp4_W.
get( 2 ), -vp4_W.
get( 3 ) );
151 cosL = vp4_Lepp.
dot( vp4_W ) / ( vp4_Lepp.
d3mag() * vp4_W.
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;
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 );
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 );
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;
189 for (
int index = 0; index < nAmps; index++ )
191 switch ( type[index] )
195 NRS( m,
q, rS, rS1, a_delta, b_delta, mA, m0_1430_S, width0_1430_S, amplitude_temp,
197 amplitude[index] = amplitude_temp;
198 delta[index] = delta_temp;
204 ResonanceP( m,
q, mV, mA, V_0, A1_0, A2_0, m0, width0, rBW, amplitude_temp, delta_temp,
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;
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;
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;
240 std::cout <<
"No such form factor type!!!" << std::endl;
247 double I, I1, I2, I3, I4, I5, I6, I7, I8, I9;
249 double cosV2 = cosV * cosV;
250 double sinV = sqrt( 1.0 - cosV2 );
251 double sinV2 = sinV * sinV;
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;
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 );
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;
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;
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;
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,
286 double pKPi = getPStar( mD, m,
q );
287 double mD2 = mD * mD;
289 double m02 = m0 * m0;
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;
301 double H0 = 0.5 / ( m *
q ) *
302 ( ( mD2 -
m2 - q2 ) * summDm * A1 - 4.0 * ( mD2 * pKPi * pKPi ) / summDm * A2 );
307 double B_Kstar = 2. / 3.;
308 double pStar0 = getPStar( m0, mPi, mK );
309 double alpha = sqrt( 3. * Pi * B_Kstar / ( pStar0 * width0 ) );
312 double F = getF1( m, m0, mPi, mK, rBW );
313 double width = getWidth1( m, m0, mPi, mK, width0, rBW );
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 ) );
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 );
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,
331 static const double tmp = ( mK + mPi ) * ( mK + mPi );
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 );
346 x = sqrt(
m2 / tmp - 1.0 );
351 x = sqrt( m2_K0_1430 / tmp - 1.0 );
353 Pm *= m_K0_1430 * width_K0_1430 /
354 sqrt( ( m2_K0_1430 -
m2 ) * ( m2_K0_1430 -
m2 ) + m2_K0_1430 * width * width );
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 );
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;
366 EvtComplex ci(
cos( delta ),
sin( delta ) );
367 EvtComplex amp = ci * rS * Pm;
370 F10 = amp * pKPi * mD / ( 1. - q2 / mA2 );
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,
377 double pKPi = getPStar( mD, m,
q );
378 double mD2 = mD * mD;
380 double m02 = m0 * m0;
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 ) );
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;
396 EvtComplex amp =
C / EvtComplex( AA, BB );
397 amplitude =
abs( amp );
398 delta = atan2(
imag( amp ),
real( amp ) );
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;
406double EvtDToKpienu::getPStar(
double m,
double m1,
double m2 ) {
410 double x =
s + s1 - s2;
411 double t = 0.25 *
x *
x /
s - s1;
413 if (
t > 0.0 ) { p = sqrt(
t ); }
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;
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 );
442 1. / sqrt( ( rBW2 * pStar02 - 3. ) * ( rBW2 * pStar02 - 3. ) + 9. * rBW2 * pStar02 );
443 double F = pStar2 / pStar02 *
B / B0;
447double EvtDToKpienu::getWidth0(
double m,
double m0,
double m_c1,
double m_c2,
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;
455double EvtDToKpienu::getWidth1(
double m,
double m0,
double m_c1,
double m_c2,
double width0,
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;
464double EvtDToKpienu::getWidth2(
double m,
double m0,
double m_c1,
double m_c2,
double width0,
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;
473EvtComplex EvtDToKpienu::getCoef(
double rho,
double phi ) {
474 EvtComplex coef( rho *
cos( phi ), rho *
sin( phi ) );
Evt3Rank3C conj(const Evt3Rank3C &t2)
double imag(const EvtComplex &c)
double abs2(const EvtComplex &c)
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
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
void getName(std::string &name)
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)
void checkNArg(int a1, int a2=-1, int a3=-1, int a4=-1)
void setProb(double prob)
static int getStdHep(EvtId id)
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
EvtVector4R cross(const EvtVector4R &v2)
void set(int i, double d)