Line data Source code
1 : // SusyWidthFunctions.cc is a part of the PYTHIA event generator.
2 : // Copyright (C) 2015 Torbjorn Sjostrand
3 : // Authors: N. Desai, P. Skands
4 : // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
5 : // Please respect the MCnet Guidelines, see GUIDELINES for details.
6 :
7 : // Function definitions (not found in the header) for the
8 : // SUSY Resonance three-body decay width classes.
9 :
10 : #include "Pythia8/SusyResonanceWidths.h"
11 : #include "Pythia8/SusyWidthFunctions.h"
12 : #include "Pythia8/SusyCouplings.h"
13 : #include "Pythia8/PythiaComplex.h"
14 :
15 : namespace Pythia8 {
16 :
17 : //==========================================================================
18 :
19 : // The WidthFunctions class.
20 : // Functions to be integrated for calculating the 3-body decay widths.
21 :
22 : //--------------------------------------------------------------------------
23 :
24 : void WidthFunction::setPointers( ParticleData* particleDataPtrIn,
25 : CoupSUSY* coupSUSYPtrIn, Info* infoPtrIn) {
26 :
27 0 : particleDataPtr = particleDataPtrIn;
28 0 : coupSUSYPtr = coupSUSYPtrIn;
29 0 : infoPtr = infoPtrIn;
30 0 : }
31 :
32 : //--------------------------------------------------------------------------
33 :
34 : double WidthFunction::function(double) {
35 :
36 0 : cout << "Warning using dummy width function" << endl;
37 0 : return 0.;
38 : }
39 :
40 : //--------------------------------------------------------------------------
41 :
42 : // Adapted from the CERNLIB DGAUSS routine by K.S. Kolbig.
43 :
44 : double WidthFunction::integrateGauss(double xlo, double xhi, double tol) {
45 :
46 : // 8-point unweighted.
47 : static double x8[4]={0.96028985649753623,
48 : 0.79666647741362674,
49 : 0.52553240991632899,
50 : 0.18343464249564980};
51 : static double w8[4]={0.10122853629037626,
52 : 0.22238103445337447,
53 : 0.31370664587788729,
54 : 0.36268378337836198};
55 : // 16-point unweighted.
56 : static double x16[8]={0.98940093499164993,
57 : 0.94457502307323258,
58 : 0.86563120238783174,
59 : 0.75540440835500303,
60 : 0.61787624440264375,
61 : 0.45801677765722739,
62 : 0.28160355077925891,
63 : 0.09501250983763744};
64 : static double w16[8]={0.027152459411754095,
65 : 0.062253523938647893,
66 : 0.095158511682492785,
67 : 0.12462897125553387,
68 : 0.14959598881657673,
69 : 0.16915651939500254,
70 : 0.18260341504492359,
71 : 0.18945061045506850};
72 : // Boundary checks.
73 0 : if (xlo == xhi) {
74 0 : cerr<<"xlo = xhi"<<endl;
75 0 : return 0.0;
76 : }
77 0 : if ( xlo > xhi ) {
78 0 : cerr<<" (integrateGauss:) -> xhi < xlo"<<endl;
79 0 : return 0.0;
80 : }
81 : // Initialize.
82 : double sum = 0.0;
83 0 : double c = 0.001/abs(xhi-xlo);
84 : double zlo = xlo;
85 : double zhi = xhi;
86 :
87 : bool nextbin = true;
88 :
89 0 : while ( nextbin ) {
90 :
91 0 : double zmi = 0.5*(zhi+zlo); // midpoint
92 0 : double zmr = 0.5*(zhi-zlo); // midpoint, relative to zlo
93 :
94 : // Calculate 8-point and 16-point quadratures.
95 : double s8=0.0;
96 0 : for (int i=0;i<4;i++) {
97 0 : double dz = zmr * x8[i];
98 0 : s8 += w8[i]*(function(zmi+dz) + function(zmi-dz));
99 : }
100 0 : s8 *= zmr;
101 : double s16=0.0;
102 0 : for (int i=0;i<8;i++) {
103 0 : double dz = zmr * x16[i];
104 0 : s16 += w16[i]*(function(zmi+dz) + function(zmi-dz));
105 : }
106 0 : s16 *= zmr;
107 0 : if (abs(s16-s8) < tol*(1+abs(s16))) {
108 : // Precision in this bin OK, add to cumulative and go to next.
109 : nextbin=true;
110 0 : sum += s16;
111 : // Next bin: LO = end of current, HI = end of integration region.
112 : zlo=zhi;
113 : zhi=xhi;
114 0 : if ( zlo == zhi ) nextbin = false;
115 : } else {
116 : // Precision in this bin not OK, subdivide.
117 0 : if (1.0 + c*abs(zmr) == 1.0) {
118 0 : cerr << " (integrateGauss:) too high accuracy required"<<endl;
119 : sum = 0.0 ;
120 0 : break;
121 : }
122 : zhi=zmi;
123 : nextbin=true;
124 : }
125 0 : }
126 : return sum;
127 0 : }
128 :
129 : //==========================================================================
130 :
131 : // The StauWidths class.
132 : // Width functions for 3-body stau decays.
133 :
134 : //--------------------------------------------------------------------------
135 :
136 : double StauWidths::getWidth(int idResIn, int idIn){
137 :
138 0 : setChannel(idResIn, idIn);
139 :
140 : // Calculate integration limits and return integrated width.
141 0 : if (idResIn % 2 == 0) return 0.0;
142 0 : double width = integrateGauss(0.0,1.0,1.0e-3);
143 : return width;
144 :
145 0 : }
146 :
147 : //--------------------------------------------------------------------------
148 :
149 : void StauWidths::setChannel(int idResIn, int idIn) {
150 :
151 : // Common masses.
152 0 : idRes = abs(idResIn);
153 0 : idIn = abs(idIn);
154 0 : mRes = particleDataPtr->m0(idRes);
155 0 : m1 = particleDataPtr->m0(1000022);
156 0 : m2 = particleDataPtr->m0(idIn);
157 :
158 0 : mInt = particleDataPtr->m0(15);
159 0 : gammaInt = particleDataPtr->mWidth(15);
160 :
161 : // Couplings etc.
162 0 : f0 = 92.4; //MeV
163 0 : gf = coupSUSYPtr->GF();
164 0 : delm = mRes - m1;
165 0 : cons = pow2(f0)*pow2(gf)*(pow2(delm) - pow2(m2))
166 0 : * coupSUSYPtr->V2CKMid(1,1) / (128.0 * pow(mRes*M_PI,3));
167 :
168 0 : if (idIn == 900111) wparam = 1.16;
169 0 : else if (idIn == 113) wparam = 0.808;
170 0 : else wparam = 1.0;
171 :
172 0 : double g = coupSUSYPtr->alphaEM(mRes * mRes);
173 : int ksusy = 1000000;
174 0 : int isl = (abs(idRes)/ksusy == 2) ? (abs(idRes)%10+1)/2 + 3
175 : : (abs(idRes)%10+1)/2;
176 :
177 0 : gL = g * coupSUSYPtr->LsllX[isl][3][1] / ( sqrt(2.0) * coupSUSYPtr->cosb);
178 0 : gR = g * coupSUSYPtr->RsllX[isl][3][1] / ( sqrt(2.0) * coupSUSYPtr->cosb);
179 :
180 : // Set function switch and internal propagators depending on decay product.
181 0 : if (idIn == 111) fnSwitch = 1;
182 0 : else if (idIn == 900111 || idIn == 113) fnSwitch = 2;
183 0 : else if (idIn == 14 || idIn == 12) {
184 0 : m2 = particleDataPtr->m0(idIn-1);
185 0 : fnSwitch = 3;
186 0 : }
187 : else {
188 0 : stringstream mess;
189 0 : mess << " unknown decay channel idIn = " << idIn;
190 0 : infoPtr->errorMsg("Warning in StauWidths::setChannel:", mess.str() );
191 0 : }
192 :
193 : return;
194 0 : }
195 :
196 : //------------------------------------------------------------------------
197 :
198 : double StauWidths::function(double x){
199 :
200 : // Decay width functions documented in arXiv:1212.2886 Citron et. al.
201 : double value = 0.0;
202 0 : double qf2 = pow2(delm) - (pow2(delm) - pow2(m2)) * x;
203 0 : double fac = 1.0 / pow3(mRes);
204 0 : double term3 = (norm(gL) * qf2 + norm(gR) * mInt * mInt)
205 0 : * (delm * delm + 2.0 * m1 * delm - qf2);
206 0 : double term4 = -2.0 * real(gL * conj(gR)) * m2 * mInt * qf2;
207 :
208 : // ~tau -> pi0 nu_tau ~chi0.
209 0 : if (fnSwitch == 1 ) {
210 0 : fac *= pow2(delm) - pow2(m2);
211 0 : double term1 = sqrt((pow2(delm) - qf2) * (pow2(delm + 2 * m1) - qf2));
212 0 : double term2 = pow2(qf2 - pow2(m2)) / qf2 / (pow2(qf2 - pow2(mInt))
213 0 : + pow2(mInt*gammaInt));
214 0 : value = fac * (term1 * term2 * (term3 + term4));
215 0 : }
216 :
217 0 : else if (fnSwitch == 2 ) {
218 0 : double term1 = sqrt((pow2(delm) - qf2) * (pow2(delm + 2 * m1) - qf2));
219 0 : double term2 = pow2(qf2 - pow2(m2)) * (qf2 + pow2(m2))
220 0 : / (qf2 * qf2 * (pow2(qf2 - pow2(mInt)) + pow2(mInt*gammaInt)));
221 0 : value = fac * (term1 * term2 * (term3 + term4));
222 0 : }
223 :
224 0 : else if (fnSwitch == 3 ) {
225 0 : double qf4 = qf2 * qf2;
226 0 : double m24 = pow2(m2*m2);
227 0 : double term1 = sqrt((pow2(delm) - qf2) * (pow2(delm + 2 * m1) - qf2));
228 0 : double term2a = 1.0 / (pow2(qf2 - pow2(mInt)) + pow2(mInt*gammaInt)) / qf4;
229 0 : double term2b = 12 * m24 * qf4 * log(qf2/pow2(m2))
230 0 : + (qf4 - m24) * (qf4 - 8 * m2 * m2 * qf2 + m24);
231 0 : value = fac * (term1 * term2a * term2b * (term3 + term4));
232 0 : }
233 :
234 : else {
235 0 : stringstream mess;
236 0 : mess << " unknown decay channel fnSwitch = " << fnSwitch;
237 0 : infoPtr->errorMsg("Warning in StauWidths::function:", mess.str() );
238 0 : }
239 :
240 0 : return value;
241 :
242 0 : }
243 :
244 : //==========================================================================
245 :
246 : } //end namespace Pythia8
|
