20#include "EvtGenBase/EvtPatches.hh"
21#include "EvtGenBase/EvtParticle.hh"
22#include "EvtGenBase/EvtGenKine.hh"
23#include "EvtGenBase/EvtPDL.hh"
24#include "EvtGenBase/EvtReport.hh"
25#include "EvtGenModels/EvtDTopipienu.hh"
26#include "EvtGenBase/EvtComplex.hh"
27#include "EvtGenBase/EvtDecayTable.hh"
33 model_name=
"DTopipienu";
52 std::cout <<
"EvtDTopipienu ==> Initialization !" << std::endl;
55 if ( parnum == D0 || parnum == D0B ) {
59 }
else if (parnum == DP || parnum == DM ) {
84 width0_omega = 0.00849;
100 Pi = atan2(0.0,-1.0);
102 root2d3 = sqrt(2./3);
112 cout <<
"EvtDTopipienu: setProbMax = " << ProbMax << endl;
148 if(pid == -211) charm = 1;
150 double m2, q2, cosV, cosL, chi;
151 KinVGen(pi1, pi2, e, nu, charm, m2, q2, cosV, cosL, chi);
152 double prob = calPDF(m2, q2, cosV, cosL, chi);
161 m2 = vp4_KPi.
mass2();
165 boost.
set(vp4_KPi.
get(0), -vp4_KPi.
get(1), -vp4_KPi.
get(2), -vp4_KPi.
get(3));
178 double sinx =
C.cross(V).dot(D);
179 double cosx =
C.dot(D);
180 chi = sinx>0? acos(cosx): -acos(cosx);
181 if(charm==-1) chi=-chi;
184double EvtDTopipienu::calPDF(
double m2,
double q2,
double cosV,
double cosL,
double chi) {
204 for(
int index=first; index<last; index++)
210 ResonanceGS(m,
q, D0_mD, D0_mPi1, D0_mPi2, f11, f21, f31);
211 coef = getCoef(rho, phi);
219 ResonanceGS(m,
q, Dp_mD, Dp_mPi1, Dp_mPi2, f11, f21, f31);
220 coef = getCoef(rho, phi);
228 ResonancePGScbw(m,
q, f11, f21, f31);
229 coef = getCoef(rho_omega, phi_omega);
237 ResonanceSBugg(m,
q, f10);
238 coef = getCoef(rho_S, phi_S);
239 F10 = F10 + coef*f10;
244 std::cout<<
"No such form factor type!!!"<<std::endl;
251 double I,I1,I2,I3,I4,I5,I6,I7,I8,I9;
253 double cosV2 = cosV*cosV;
254 double sinV = sqrt(1.0-cosV2);
255 double sinV2 = sinV*sinV;
257 EvtComplex F1 = F10 + F11*cosV + F12*(1.5*cosV2-0.5);
258 EvtComplex F2 = F21*root1d2 + F22*cosV*root3d2;
259 EvtComplex F3 = F31*root1d2 + F32*cosV*root3d2;
261 I1 = 0.25*(
abs2(F1) + 1.5*sinV2*(
abs2(F2) +
abs2(F3) ) );
262 I2 = -0.25*(
abs2(F1) - 0.5*sinV2*(
abs2(F2) +
abs2(F3) ) );
263 I3 = -0.25*(
abs2(F2) -
abs2(F3) )*sinV2;
269 I9 =
imag(
conj(F3)*F2 )*sinV2*(-0.5);
271 double coschi =
cos(chi);
272 double sinchi =
sin(chi);
273 double sin2chi = 2.0*sinchi*coschi;
274 double cos2chi = 1.0 - 2.0*sinchi*sinchi;
276 double sinL = sqrt(1.-cosL*cosL);
277 double sinL2 = sinL*sinL;
278 double sin2L = 2.0*sinL*cosL;
279 double cos2L = 1.0 - 2.0*sinL2;
281 I = I1 + I2*cos2L + I3*sinL2*cos2chi + I4*sin2L*coschi + I5*sinL*coschi + I6*cosL + I7*sinL*sinchi + I8*sin2L*sinchi + I9*sinL2*sin2chi;
287 double pKPi = getPStar(massD, m,
q);
288 double mD2 = massD*massD;
294 double summDm = massD+m;
295 double V = V_0 /(1.0-q2/(mV2));
296 double A1 = A1_0/(1.0-q2/(mA2));
297 double A2 = A2_0/(1.0-q2/(mA2));
298 double A = summDm*A1;
299 double B = 2.0*massD*pKPi/summDm*V;
302 double H0 = 0.5/(m*
q)*((mD2-m2-q2)*summDm*A1-4.0*(mD2*pKPi*pKPi)/summDm*A2);
308 double pStar = getPStar(m, massPi1, massPi2);
309 double pStar0 = getPStar(m0, massPi1, massPi2);
310 double alpha = sqrt(3.0*Pi*BF/(pStar0*width0));
313 double F = getF1(m, m0, massPi1, massPi2, rBW);
314 EvtComplex C(m02*(1.0+width0*getGx(m0, pStar0, massPi1, massPi2)/m0)*F, 0.0);
315 double AA = m02-m2+width0*getFx(m02, m2, pStar, pStar0, massPi1, massPi2);
316 double BB = -m0*getWidthrho(m, m0, width0, pStar, pStar0);
320 double alpham2 =
alpha*2.0;
321 F11 = amp*alpham2*
q*H0*root2;
322 F21 = amp*alpham2*
q*(Hp+Hm);
323 F31 = amp*alpham2*
q*(Hp-Hm);
328 double pKPi = getPStar(Dp_mD, m,
q);
329 double mD2 = Dp_mD*Dp_mD;
331 double m02 = m0_omega*m0_omega;
336 double summDm = Dp_mD+m;
337 double V = V_0 /(1.0-q2/(mV2));
338 double A1 = A1_0/(1.0-q2/(mA2));
339 double A2 = A2_0/(1.0-q2/(mA2));
340 double A = summDm*A1;
341 double B = 2.0*Dp_mD*pKPi/summDm*V;
343 double H0 = 0.5/(m*
q)*((mD2-m2-q2)*summDm*A1-4.0*(mD2*pKPi*pKPi)/summDm*A2);
349 double pStar = getPStar(m, Dp_mPi1, Dp_mPi2);
350 double pStar0 = getPStar(m0_omega, Dp_mPi1, Dp_mPi2);
351 double alpha = sqrt(3.0*Pi*BF_omega/(pStar0*width0_omega));
354 double F = getF1(m, m0_omega, Dp_mPi1, Dp_mPi2, rBW);
360 double BB1 = -m0_omega*width0_omega;
364 pStar0 = getPStar(m0, Dp_mPi1, Dp_mPi2);
365 EvtComplex C2(mR2*(1.0+width0*getGx(m0, pStar0, Dp_mPi1, Dp_mPi2)/m0), 0.0);
366 double AA2 = mR2-m2+width0*getFx(mR2, m2, pStar, pStar0, Dp_mPi1, Dp_mPi2);
367 double BB2 = -m0*getWidthrho(m, m0, width0, pStar, pStar0);
372 double alpham2 =
alpha*2.0;
373 F11 = amp*alpham2*
q*H0*root2;
374 F21 = amp*alpham2*
q*(Hp+Hm);
375 F31 = amp*alpham2*
q*(Hp-Hm);
378void EvtDTopipienu::ResonanceSBugg(
double m,
double q,
EvtComplex& F10)
380 double pKPi = getPStar(Dp_mD, m,
q);
385 double sA = 0.41*mPi*mPi;
387 double mr2 = m0_S*m0_S;
398 Gamma1 = getrho(m2,mPi)*getG1(m2,mr)*(m2-sA)/(mr2-sA);
399 Gamma2 = getrho(m2,mKa)*0.6*getG1(m2,mr)*(m2/mr2)*
exp(-
alpha*fabs(m2-4.0*mKa*mKa));
400 Gamma3 = getrho(m2,mEt)*0.2*getG1(m2,mr)*(m2/mr2)*
exp(-
alpha*fabs(m2-4.0*mEt*mEt));
401 if (m > 4*mPi) Gamma4 = ciR*mr*g4pi*(1.0+
exp(7.082-2.845*mr2))/(1.0+
exp(7.082-2.845*m2));
403 double AA = mr2-m2-getG1(m2,mr)*(m2-sA)*getZ(m2, mr2)/(mr2-sA);
404 EvtComplex amp = ciR/(ciR*AA-ciM*(Gamma1+Gamma2+Gamma3+Gamma4));
405 F10 = amp*pKPi*Dp_mD/(1.0-q2/mA2);
408double EvtDTopipienu::getPStar(
double m,
double m1,
double m2)
413 double x =
s + s1 - s2;
414 double t = 0.25*
x*
x/
s - s1;
419 std::cout <<
" Hello, pstar is less than 0.0" << std::endl;
425double EvtDTopipienu::getF1(
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(1.+rBW2*pStar2);
433 double B0 = 1./sqrt(1.+rBW2*pStar02);
434 double F = pStar/pStar0*
B/B0;
438double EvtDTopipienu::getF2(
double m,
double m0,
double m_c1,
double m_c2,
double rBW)
440 double pStar = getPStar(m, m_c1, m_c2);
441 double pStar0 = getPStar(m0, m_c1, m_c2);
442 double rBW2 = rBW*rBW;
443 double pStar2 = pStar*pStar;
444 double pStar02 = pStar0*pStar0;
445 double B = 1./sqrt((rBW2*pStar2-3.)*(rBW2*pStar2-3.)+9.*rBW2*pStar2);
446 double B0 = 1./sqrt((rBW2*pStar02-3.)*(rBW2*pStar02-3.)+9.*rBW2*pStar02);
447 double F = pStar2/pStar02*
B/B0;
451double EvtDTopipienu::getWidth0(
double m,
double m0,
double m_c1,
double m_c2,
double width0)
453 double pStar = getPStar(m, m_c1, m_c2);
454 double pStar0 = getPStar(m0, m_c1, m_c2);
455 double width = width0*pStar/pStar0*m0/m;
459double EvtDTopipienu::getWidth1(
double m,
double m0,
double m_c1,
double m_c2,
double width0,
double rBW)
461 double pStar = getPStar(m, m_c1, m_c2);
462 double pStar0 = getPStar(m0, m_c1, m_c2);
463 double F = getF1(m, m0, m_c1, m_c2, rBW);
464 double width = width0*pStar/pStar0*m0/m*F*F;
468double EvtDTopipienu::getWidth2(
double m,
double m0,
double m_c1,
double m_c2,
double width0,
double rBW)
470 double pStar = getPStar(m, m_c1, m_c2);
471 double pStar0 = getPStar(m0, m_c1, m_c2);
472 double F = getF2(m, m0, m_c1, m_c2, rBW);
473 double width = width0*pStar/pStar0*m0/m*F*F;
477EvtComplex EvtDTopipienu::getCoef(
double rho,
double phi)
483inline double EvtDTopipienu::getGx(
double m0,
double p0,
double m_c1,
double m_c2)
486 double MPI = 0.5*(m_c1+m_c2);
487 Gg = 3*MPI*MPI*log((m0+2*p0)/(2*MPI))/(Pi*p0*p0) + m0/(2*Pi*p0) - MPI*MPI*m0/(Pi*p0*p0*p0);
491inline double EvtDTopipienu::getFx(
double mr2,
double sx,
double p,
double p0,
double m_c1,
double m_c2)
494 Fx = mr2/(pow(p0,3))*(p*p*(getHx(sx,p,m_c1,m_c2)-getHx(mr2,p0,m_c1,m_c2))+(mr2-sx)*p0*p0*getdh(mr2,p0,m_c1,m_c2));
498inline double EvtDTopipienu::getHx(
double sx,
double p,
double m_c1,
double m_c2)
502 Hx = 2*p*log((m+2*p)/(m_c1+m_c2))/(Pi*m);
506inline double EvtDTopipienu::getdh(
double mr2,
double p0,
double m_c1,
double m_c2)
508 double mass = sqrt(mr2);
513inline double EvtDTopipienu::getG1(
double m2,
double Mr)
519 double gg1 = Mr*(b1+b2*m2)*
exp(-(m2-Mr2)/
A);
523inline double EvtDTopipienu::getZ(
double m2,
double Mr2)
525 double zz = (getRho(m2,mPi)*log((1.0-getRho(m2,mPi))/(1.0+getRho(m2,mPi)))
526 -getRho(Mr2,mPi)*log((1.0-getRho(Mr2,mPi))/(1.0+getRho(Mr2,mPi))))/Pi;
531inline double EvtDTopipienu::getRho(
double m2,
double mX)
534 if ((1.0-4.0*mX*mX/m2)>0)
535 rho = sqrt(1.0-4.0*mX*mX/m2);
540inline EvtComplex EvtDTopipienu::getrho(
double m2,
double mX)
545 if ((1.0-4.0*mX*mX/m2)>0)
546 rho = ciR*sqrt(1.0-4.0*mX*mX/m2);
547 else rho = ciM*sqrt(4.0*mX*mX/m2-1.0);
550inline double EvtDTopipienu::getWidthrho(
double m,
double m0,
double width0,
double p,
double p0 )
552 double widthRho = 0.0;
553 widthRho = width0*pow(p/p0, 3)*(m0/m);
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
Evt3Rank3C conj(const Evt3Rank3C &t2)
double imag(const EvtComplex &c)
double abs2(const EvtComplex &c)
EvtComplex exp(const EvtComplex &c)
EvtDiracSpinor boostTo(const EvtDiracSpinor &sp, const EvtVector4R p4)
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)
static EvtId getId(const 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)
double dot(const EvtVector4R &v2) const
EvtVector4R cross(const EvtVector4R &v2)
void set(int i, double d)