33 model_name=
"DToKpienu";
47 std::cout <<
"EvtDToKpienu ==> Initialization !" << std::endl;
55 width0_1430_S = 0.270;
131 if(pid == -321) charm = 1;
133 double m2, q2, cosV, cosL, chi;
134 KinVGen(K,
pi, e, nu, charm, m2, q2, cosV, cosL, chi);
135 double prob = calPDF(m2, q2, cosV, cosL, chi);
145 m2 = vp4_KPi.
mass2();
149 boost.
set(vp4_KPi.
get(0), -vp4_KPi.
get(1), -vp4_KPi.
get(2), -vp4_KPi.
get(3));
162 double sinx =
C.cross(V).dot(D);
163 double cosx =
C.dot(D);
164 chi = sinx>0? acos(cosx): -acos(cosx);
165 if(charm==-1) chi=-chi;
168double EvtDToKpienu::calPDF(
double m2,
double q2,
double cosV,
double cosL,
double chi) {
169 double delta[5] = {0};
170 double amplitude[5] = {0};
191 double amplitude_temp, delta_temp;
193 for(
int index=0; index<nAmps; index++)
199 NRS(m,
q, rS, rS1, a_delta, b_delta, mA, m0_1430_S, width0_1430_S, amplitude_temp, delta_temp, f10);
200 amplitude[index] = amplitude_temp;
201 delta[index] = delta_temp;
207 ResonanceP(m,
q, mV, mA, V_0, A1_0, A2_0, m0, width0, rBW, amplitude_temp, delta_temp, f11, f21, f31);
208 amplitude[index] = amplitude_temp;
209 delta[index] = delta_temp;
210 coef = getCoef(rho, phi);
219 ResonanceP(m,
q, mV, mA, V_0, A1_0, A2_0, m0_1410, width0_1410, rBW, amplitude_temp, delta_temp, f11, f21, f31);
220 amplitude[index] = amplitude_temp;
221 delta[index] = delta_temp;
222 coef = getCoef(rho_1410, phi_1410);
230 ResonanceD(m,
q, mV, mA, TV_0, T1_0, T2_0, m0_1430, width0_1430, rBW, amplitude_temp, delta_temp, f12, f22, f32);
231 amplitude[index] = amplitude_temp;
232 delta[index] = delta_temp;
233 coef = getCoef(rho_1430, phi_1430);
241 std::cout<<
"No such form factor type!!!"<<std::endl;
248 double I,I1,I2,I3,I4,I5,I6,I7,I8,I9;
250 double cosV2 = cosV*cosV;
251 double sinV = sqrt(1.0-cosV2);
252 double sinV2 = sinV*sinV;
254 EvtComplex F1 = F10 + F11*cosV + F12*(1.5*cosV2-0.5);
255 EvtComplex F2 = F21*root1d2 + F22*cosV*root3d2;
256 EvtComplex F3 = F31*root1d2 + F32*cosV*root3d2;
258 I1 = 0.25*(
abs2(F1) + 1.5*sinV2*(
abs2(F2) +
abs2(F3) ) );
259 I2 = -0.25*(
abs2(F1) - 0.5*sinV2*(
abs2(F2) +
abs2(F3) ) );
260 I3 = -0.25*(
abs2(F2) -
abs2(F3) )*sinV2;
266 I9 =
imag(
conj(F3)*F2 )*sinV2*(-0.5);
268 double coschi =
cos(chi);
269 double sinchi =
sin(chi);
270 double sin2chi = 2.0*sinchi*coschi;
271 double cos2chi = 1.0 - 2.0*sinchi*sinchi;
273 double sinL = sqrt(1.-cosL*cosL);
274 double sinL2 = sinL*sinL;
275 double sin2L = 2.0*sinL*cosL;
276 double cos2L = 1.0 - 2.0*sinL2;
278 I = I1 + I2*cos2L + I3*sinL2*cos2chi + I4*sin2L*coschi + I5*sinL*coschi + 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,
double A1_0,
double A2_0,
double m0,
double width0,
double rBW,
double& amplitude,
double&
delta,
EvtComplex& F11,
EvtComplex& F21,
EvtComplex& F31)
284 double pKPi = getPStar(mD, m,
q);
291 double summDm = mD+m;
292 double V = V_0 /(1.0-q2/(mV2));
293 double A1 = A1_0/(1.0-q2/(mA2));
294 double A2 = A2_0/(1.0-q2/(mA2));
295 double A = summDm*A1;
296 double B = 2.0*mD*pKPi/summDm*V;
299 double H0 = 0.5/(m*
q)*((mD2-m2-q2)*summDm*A1-4.0*(mD2*pKPi*pKPi)/summDm*A2);
304 double B_Kstar = 2./3.;
305 double pStar0 = getPStar(m0, mPi, mK);
306 double alpha = sqrt(3.*Pi*B_Kstar/(pStar0*width0));
309 double F = getF1(m, m0, mPi, mK, rBW);
310 double width = getWidth1(m, m0, mPi, mK, width0, rBW);
314 double BB = -m0*width;
316 amplitude =
abs(amp);
319 double alpham2 =
alpha*2.0;
320 F11 = amp*alpham2*
q*H0*root2;
321 F21 = amp*alpham2*
q*(Hp+Hm);
322 F31 = amp*alpham2*
q*(Hp-Hm);
325void EvtDToKpienu::NRS(
double m,
double q,
double rS,
double rS1,
double a_delta,
double b_delta,
double mA,
double m0,
double width0,
double& amplitude,
double&
delta,
EvtComplex& F10)
327 static const double tmp = (mK+mPi)*(mK+mPi);
332 double pKPi = getPStar(mD, m,
q);
333 double m_K0_1430 = m0;
334 double width_K0_1430 = width0;
335 double m2_K0_1430 = m_K0_1430*m_K0_1430;
336 double width = getWidth0(m, m_K0_1430, mPi, mK, width_K0_1430);
341 x = sqrt(m2/tmp-1.0);
344 x = sqrt(m2_K0_1430/tmp-1.0);
346 Pm *= m_K0_1430*width_K0_1430/sqrt((m2_K0_1430-m2)*(m2_K0_1430-m2)+m2_K0_1430*width*width);
350 double pStar = getPStar(m, mPi, mK);
351 double delta_bg = atan(2.*a_delta*pStar/(2.+a_delta*b_delta*pStar*pStar));
352 delta_bg = (delta_bg>0)? delta_bg: (delta_bg+Pi);
354 double delta_K0_1430 = atan(m_K0_1430*width/(m2_K0_1430-m2));
355 delta_K0_1430 = (delta_K0_1430>0)? delta_K0_1430: (delta_K0_1430+Pi);
356 delta = delta_bg + delta_K0_1430;
362 F10 = amp*pKPi*mD/(1.-q2/mA2);
365void EvtDToKpienu::ResonanceD(
double m,
double q,
double mV,
double mA,
double TV_0,
double T1_0,
double T2_0,
double m0,
double width0,
double rBW,
double& amplitude,
double&
delta,
EvtComplex& F12,
EvtComplex& F22,
EvtComplex& F32)
367 double pKPi = getPStar(mD, m,
q);
374 double summDm = mD+m;
375 double TV = TV_0 /(1.0-q2/(mV2));
376 double T1 = T1_0/(1.0-q2/(mA2));
377 double T2 = T2_0/(1.0-q2/(mA2));
380 double F = getF2(m, m0, mPi, mK, rBW);
381 double width = getWidth2(m, m0, mPi, mK, width0, rBW);
384 double BB = -m0*width;
387 amplitude =
abs(amp);
390 F12 = amp*mD*pKPi/3.*((mD2-m2-q2)*summDm*T1-mD2*pKPi*pKPi/summDm*T2);
391 F22 = amp*root2d3*mD*m*
q*pKPi*summDm*T1;
392 F32 = amp*root2d3*2.*mD2*m*
q*pKPi*pKPi/summDm*TV;
395double EvtDToKpienu::getPStar(
double m,
double m1,
double m2)
400 double x =
s + s1 - s2;
401 double t = 0.25*
x*
x/
s - s1;
406 std::cout <<
" Hello, pstar is less than 0.0" << std::endl;
412double EvtDToKpienu::getF1(
double m,
double m0,
double m_c1,
double m_c2,
double rBW)
414 double pStar = getPStar(m, m_c1, m_c2);
415 double pStar0 = getPStar(m0, m_c1, m_c2);
416 double rBW2 = rBW*rBW;
417 double pStar2 = pStar*pStar;
418 double pStar02 = pStar0*pStar0;
419 double B = 1./sqrt(1.+rBW2*pStar2);
420 double B0 = 1./sqrt(1.+rBW2*pStar02);
421 double F = pStar/pStar0*
B/B0;
425double EvtDToKpienu::getF2(
double m,
double m0,
double m_c1,
double m_c2,
double rBW)
427 double pStar = getPStar(m, m_c1, m_c2);
428 double pStar0 = getPStar(m0, m_c1, m_c2);
429 double rBW2 = rBW*rBW;
430 double pStar2 = pStar*pStar;
431 double pStar02 = pStar0*pStar0;
432 double B = 1./sqrt((rBW2*pStar2-3.)*(rBW2*pStar2-3.)+9.*rBW2*pStar2);
433 double B0 = 1./sqrt((rBW2*pStar02-3.)*(rBW2*pStar02-3.)+9.*rBW2*pStar02);
434 double F = pStar2/pStar02*
B/B0;
438double EvtDToKpienu::getWidth0(
double m,
double m0,
double m_c1,
double m_c2,
double width0)
440 double pStar = getPStar(m, m_c1, m_c2);
441 double pStar0 = getPStar(m0, m_c1, m_c2);
442 double width = width0*pStar/pStar0*m0/m;
446double EvtDToKpienu::getWidth1(
double m,
double m0,
double m_c1,
double m_c2,
double width0,
double rBW)
448 double pStar = getPStar(m, m_c1, m_c2);
449 double pStar0 = getPStar(m0, m_c1, m_c2);
450 double F = getF1(m, m0, m_c1, m_c2, rBW);
451 double width = width0*pStar/pStar0*m0/m*F*F;
455double EvtDToKpienu::getWidth2(
double m,
double m0,
double m_c1,
double m_c2,
double width0,
double rBW)
457 double pStar = getPStar(m, m_c1, m_c2);
458 double pStar0 = getPStar(m0, m_c1, m_c2);
459 double F = getF2(m, m0, m_c1, m_c2, rBW);
460 double width = width0*pStar/pStar0*m0/m*F*F;
464EvtComplex EvtDToKpienu::getCoef(
double rho,
double phi)
double sin(const BesAngle a)
double cos(const BesAngle a)
Evt3Rank3C conj(const Evt3Rank3C &t2)
double imag(const EvtComplex &c)
double abs2(const EvtComplex &c)
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
****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)