5 #include "Pythia8/PhaseSpace.h"
10 class WprimeWZFlatPtMass6;
30 const PhaseSpace* phaseSpacePtr,
33 if (sigmaProcessPtr->nFinal() != 1)
return 0.;
40 int idRes = sigmaProcessPtr->resonanceA();
41 double mRes = particleDataPtr->m0(idRes);
42 double wRes = particleDataPtr->mWidth(idRes);
43 double m2Res = mRes*mRes;
44 double gamMRat = wRes/mRes;
45 double sHat = phaseSpacePtr->sHat();
46 double weightBW = m2Res*m2Res + sHat*sHat*(1 + gamMRat*gamMRat) - 2.*sHat*m2Res;
47 double rH = sqrt(sHat);
66 bool isVetoed =
false;
72 if (
v.idAbs() != 34)
continue;
77 if(d_W.idAbs() !=24 ){
83 Vec4 pv_W_orig = d_W.p();
84 Vec4 pv_Z_orig = d_Z.p();
85 double pTW = sqrt(
pow(pv_W.px(),2) +
pow(pv_W.py(),2));
86 double mW = sqrt(
pow(pv_W.e(),2) -
pow(pv_W.px(),2) -
pow(pv_W.py(),2) -
pow(pv_W.pz(),2));
87 double mZ = sqrt(
pow(pv_Z.e(),2) -
pow(pv_Z.px(),2) -
pow(pv_Z.py(),2) -
pow(pv_Z.pz(),2));
90 if((
pow(mWflat,2) -
pow(mW,2))/
pow(pTW,2) < 1.0){
91 apT = sqrt(1 - (
pow(mWflat,2) -
pow(mW,2))/
pow(pTW,2));
97 pv_W.px(apT*pv_W.px());
98 pv_W.py(apT*pv_W.py());
103 pv_Z.px(apT*pv_Z.px());
104 pv_Z.py(apT*pv_Z.py());
106 d_Z.m(sqrt(
pow(pv_Z.e(),2) -
pow(pv_Z.px(),2) -
pow(pv_Z.py(),2) -
pow(pv_Z.pz(),2)));
107 double mZref = d_Z.m();
141 double mWt = ((mMax - mMin)/
log(mMax/mMin))*
log(mW/mMin)+mMin;
147 Vec4 pd1_W = d1_W.p();
148 Vec4 pd2_W = d2_W.p();
150 pd1_W.bstback(pv_W_orig);
151 pd2_W.bstback(pv_W_orig);
153 double m1_W_orig = sqrt(
pow(pd1_W.e(),2) -
pow(pd1_W.px(),2) -
pow(pd1_W.py(),2) -
pow(pd1_W.pz(),2));
154 double m2_W_orig = sqrt(
pow(pd2_W.e(),2) -
pow(pd2_W.px(),2) -
pow(pd2_W.py(),2) -
pow(pd2_W.pz(),2));
155 double mWchild = sqrt(
pow(pd1_W.e()+pd2_W.e(),2)-
pow(pd1_W.px()+pd2_W.px(),2)-
pow(pd1_W.py()+pd2_W.py(),2)-
pow(pd1_W.pz()+pd2_W.pz(),2));
157 double r_W_orig = 1 - (
pow(m1_W_orig,2) +
pow(m2_W_orig,2))/
pow(mW,2);
158 double r_W_flat = 1 - (
pow(m1_W_orig,2) +
pow(m2_W_orig,2))/
pow(mWflat,2);
159 double mR_W = (mWflat*r_W_flat)/(mWchild*r_W_orig);
162 pd1_W.px(mR_W*pd1_W.px());
163 pd1_W.py(mR_W*pd1_W.py());
164 pd1_W.pz(mR_W*pd1_W.pz());
165 pd1_W.e(sqrt(
pow(pd1_W.px(),2)+
pow(pd1_W.py(),2)+
pow(pd1_W.pz(),2)+
pow(m1_W_orig,2)));
167 pd2_W.px(mR_W*pd2_W.px());
168 pd2_W.py(mR_W*pd2_W.py());
169 pd2_W.pz(mR_W*pd2_W.pz());
170 pd2_W.e(sqrt(
pow(pd2_W.px(),2)+
pow(pd2_W.py(),2)+
pow(pd2_W.pz(),2)+
pow(m2_W_orig,2)));
173 pd1_W.bst(pv_W); d1_W.p(pd1_W);
174 pd2_W.bst(pv_W); d2_W.p(pd2_W);
175 double m1_W = sqrt(
pow(pd1_W.e(),2) -
pow(pd1_W.px(),2) -
pow(pd1_W.py(),2) -
pow(pd1_W.pz(),2));
176 double m2_W = sqrt(
pow(pd2_W.e(),2) -
pow(pd2_W.px(),2) -
pow(pd2_W.py(),2) -
pow(pd2_W.pz(),2));
187 if(mFrac < 0.0425)
return 1
e-12/(-1.293+1.098e+2*mFrac-2.800e+3*mFrac*mFrac+2.345e+4*mFrac*mFrac*mFrac);
188 if(mFrac < 0.073)
return 1.248e-12*(
exp(1.158+18.34*mFrac));
190 return 5.733e-10*
pow(mFrac,-3.798-0.6555*
log(mFrac))/
pow(1.427-mFrac,30.017);
195 double pe0 = 9.705/2000.;
196 double pe1 = -1.27668e-03;
198 double weightHighpT =1./(
exp(pe0+pe1*rH));
200 double p0 = 0.00405295;
201 double p1 = -1.15389e-06;
202 double p2 = -8.83305e-10;
203 double p3 = 1.02983e-12;
204 double p4 = -3.64486e-16;
205 double p5 = 6.05783e-20;
206 double p6 = -4.74988e-24;
207 double p7 = 1.40627e-28;
210 if(rH < 400.) weightFinal *= 0.5;
212 return weightHighpT * weightFinal;