]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PYTHIA8/pythia8140/src/SusyCouplings.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
coverity
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8140 / src / SusyCouplings.cxx
CommitLineData
b584e2f5 1// SusyCouplings.cc is a part of the PYTHIA event generator.
2// Copyright (C) 2010 Peter Skands, Torbjorn Sjostrand.
3// PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4// Please respect the MCnet Guidelines, see GUIDELINES for details.
5
6// Function definitions (not found in the header) for the
7// supersymmetric couplings class.
8
9#include "SusyCouplings.h"
10
11namespace Pythia8 {
12
13//==========================================================================
14
15// The CoupSUSY class.
16
17//--------------------------------------------------------------------------
18
19// Constants: could be changed here if desired, but normally should not.
20// These are of technical nature, as described for each.
21
22// Allow verbose printout for debug purposes.
23const bool CoupSUSY::DEBUG = false;
24
25//--------------------------------------------------------------------------
26
27// Initialize SM+SUSY couplings (only performed once).
28
29void CoupSUSY::init (SusyLesHouches* slhaPtrIn, Settings* settingsPtrIn,
30 ParticleData* particleDataPtrIn, CoupSM* coupSMPtrIn) {
31
32 // Save pointers..
33 slhaPtr = slhaPtrIn;
34 settingsPtr = settingsPtrIn;
35 particleDataPtr = particleDataPtrIn;
36 coupSMPtr = coupSMPtrIn;
37
38 // Is NMSSM switched on?
39 isNMSSM = (slhaPtr->modsel(3) != 1 ? false : true);
40 settingsPtr->flag("SLHA:NMSSM",isNMSSM);
41 int nNeut = (isNMSSM ? 5 : 4);
42 int nChar = 2;
43
44 // Initialize pole masses
45 mZpole = particleDataPtr->m0(23);
46 wZpole = particleDataPtr->mWidth(23);
47 mWpole = particleDataPtr->m0(24);
48 wWpole = particleDataPtr->mWidth(24);
49
50 // Running masses and weak mixing angle
51 // (default to pole values if no running available)
52 mW = mWpole;
53 mZ = mZpole;
54 sin2W = 1.0 - pow(mW/mZ,2);
55
56 if (slhaPtr->gauge.exists(1) && slhaPtr->gauge.exists(2)
57 && slhaPtr->hmix.exists(3)) {
58 double gp=slhaPtr->gauge(1);
59 double g =slhaPtr->gauge(2);
60 double v =slhaPtr->hmix(3);
61 mW = g * v / 2.0;
62 mZ = sqrt(pow(gp,2)+pow(g,2)) * v / 2.0;
63 double tan2W = pow2(gp)/pow2(g);
64 if (DEBUG) cout << " tan2W = " << tan2W << endl;
65 sin2W = pow2(gp)/(pow2(g)+pow2(gp));
66 }
67 sinW = sqrt(sin2W);
68 cosW = sqrt(1.0-sin2W);
69
70 // Tan(beta)
71 // By default, use the running one in HMIX (if not found, use MINPAR)
72 tanb = slhaPtr->hmix.exists(2) ? slhaPtr->hmix(2) : slhaPtr->minpar(3);
73 cosb = sqrt( 1.0 / (1.0 + tanb*tanb) );
74 sinb = sqrt(max(0.0,1.0-cosb*cosb));
75
76 // tmp : verbose output
77 if (DEBUG) {
78 cout << " sin2W(Q) = " << sin2W << " mW(Q) = " << mW
79 << " mZ(Q) = " << mZ << endl;
80 cout << " vev(Q) = " << slhaPtr->hmix(3) << " tanb(Q) = " << tanb
81 << endl;
82 for (int i=1;i<=3;i++) {
83 for (int j=1;j<=3;j++) {
84 cout << " VCKM [" << i << "][" << j << "] = "
85 << scientific << setw(10) << coupSMPtr->VCKMgen(i,j) << endl;
86 }
87 }
88 }
89
90 // Shorthand for squark mixing matrices
91 SusyLesHouches::matrixblock<6> Ru(slhaPtr->usqmix);
92 SusyLesHouches::matrixblock<6> Rd(slhaPtr->dsqmix);
93 SusyLesHouches::matrixblock<6> imRu(slhaPtr->imusqmix);
94 SusyLesHouches::matrixblock<6> imRd(slhaPtr->imusqmix);
95
96 // Construct ~g couplings
97 for (int i=1 ; i<=6 ; i++) {
98 for (int j=1 ; j<=3 ; j++) {
99 LsddG[i][j] = complex( Rd(i,j) , imRd(i,j));
100 RsddG[i][j] = complex(-Rd(i,j+3), -imRd(i,j+3));
101 LsuuG[i][j] = complex( Ru(i,j) , imRu(i,j));
102 RsuuG[i][j] = complex(-Ru(i,j+3), -imRu(i,j+3));
103 }
104 }
105
106 // Construct qqZ couplings
107 for (int i=1 ; i<=6 ; i++) {
108
109 // q[i] q[i] Z (def with extra factor 2 compared to [Okun])
110 LqqZ[i] = coupSMPtr->af(i) - 2.0*coupSMPtr->ef(i)*sin2W ;
111 RqqZ[i] = - 2.0*coupSMPtr->ef(i)*sin2W ;
112
113 // tmp: verbose output
114 if (DEBUG) {
115 cout << " LqqZ [" << i << "][" << i << "] = "
116 << scientific << setw(10) << LqqZ[i]
117 << " RqqZ [" << i << "][" << i << "] = "
118 << scientific << setw(10) << RqqZ[i] << endl;
119 }
120 }
121
122 // Construct ~q~qZ couplings
123 for (int i=1 ; i<=6 ; i++) {
124
125 // Squarks can have off-diagonal couplings as well
126 for (int j=1 ; j<=6 ; j++) {
127
128 // ~d[i] ~d[j] Z
129 LsdsdZ[i][j] = 0.0;
130 RsdsdZ[i][j] = 0.0;
131 for (int k=1;k<=3;k++) {
132 complex Rdik = complex(Rd(i,k), imRd(i,k) );
133 complex Rdjk = complex(Rd(j,k), imRd(j,k) );
134 complex Rdik3 = complex(Rd(i,k+3),imRd(i,k+3));
135 complex Rdjk3 = complex(Rd(j,k+3),imRd(j,k+3));
136 LsdsdZ[i][j] += LqqZ[1] * (Rdik*conj(Rdjk));
137 RsdsdZ[i][j] += RqqZ[1] * (Rdik3*conj(Rdjk3));
138 }
139
140 // ~u[i] ~u[j] Z
141 LsusuZ[i][j] = 0.0;
142 RsusuZ[i][j] = 0.0;
143 for (int k=1;k<=3;k++) {
144 complex Ruik = complex(Ru(i,k) ,imRu(i,k) );
145 complex Rujk = complex(Ru(j,k) ,imRu(j,k) );
146 complex Ruik3 = complex(Ru(i,k+3),imRu(i,k+3));
147 complex Rujk3 = complex(Ru(j,k+3),imRu(j,k+3));
148 LsusuZ[i][j] += LqqZ[2] * (Ruik*conj(Rujk));
149 RsusuZ[i][j] += RqqZ[2] * (Ruik3*conj(Rujk3));
150 }
151
152 // tmp: verbose output
153 if (DEBUG) {
154 if (max(abs(LsdsdZ[i][j]),abs(RsdsdZ[i][j])) > 1e-6) {
155 cout << " LsdsdZ[" << i << "][" << j << "] = "
156 << scientific << setw(10) << LsdsdZ[i][j]
157 << " RsdsdZ[" << i << "][" << j << "] = "
158 << scientific << setw(10) << RsdsdZ[i][j] << endl;
159 }
160 if (max(abs(LsusuZ[i][j]),abs(RsusuZ[i][j]))> 1e-6) {
161 cout << " LsusuZ[" << i << "][" << j << "] = "
162 << scientific << setw(10) << LsusuZ[i][j]
163 << " RsusuZ[" << i << "][" << j << "] = "
164 << scientific << setw(10) << RsusuZ[i][j] << endl;
165 }
166 }
167
168 }
169
170 }
171
172 // Construct udW couplings
173
174 // Loop over up [i] and down [j] quark generation
175 for (int i=1;i<=3;i++) {
176 for (int j=1;j<=3;j++) {
177
178 // CKM matrix (use Pythia one if no SLHA)
179 // (NB: could also try input one if no running one found, but
180 // would then need to compute from Wolfenstein)
181 complex Vij=coupSMPtr->VCKMgen(i,j);
182 if (slhaPtr->vckm.exists()) {
183 Vij=complex(slhaPtr->vckm(i,j),slhaPtr->imvckm(i,j));
184 }
185
186 // u[i] d[j] W
187 LudW[i][j] = sqrt(2.0) * cosW * Vij;
188 RudW[i][j] = 0.0;
189
190 // tmp: verbose output
191 if (DEBUG) {
192 cout << " LudW [" << i << "][" << j << "] = "
193 << scientific << setw(10) << LudW[i][j]
194 << " RudW [" << i << "][" << j << "] = "
195 << scientific << setw(10) << RudW[i][j] << endl;
196 }
197 }
198 }
199
200 // Construct ~u~dW couplings
201
202 // Loop over ~u[k] and ~d[l] flavours
203 for (int k=1;k<=6;k++) {
204 for (int l=1;l<=6;l++) {
205
206 LsusdW[k][l]=0.0;
207 RsusdW[k][l]=0.0;
208
209 // Loop over u[i] and d[j] flavours
210 for (int i=1;i<=3;i++) {
211 for (int j=1;j<=3;j++) {
212
213 // CKM matrix (use Pythia one if no SLHA)
214 // (NB: could also try input one if no running one found, but
215 // would then need to compute from Wolfenstein)
216 complex Vij=coupSMPtr->VCKMgen(i,j);
217 if (slhaPtr->vckm.exists()) {
218 Vij=complex(slhaPtr->vckm(i,j),slhaPtr->imvckm(i,j));
219 }
220
221 // ~u[k] ~d[l] W (add one term for each quark flavour i,j)
222 complex Ruki = complex(Ru(k,i),imRu(k,i));
223 complex Rdlj = complex(Rd(l,j),imRd(l,j));
224 LsusdW[k][l] += sqrt(2.0) * cosW * Vij * Ruki * conj(Rdlj);
225 RsusdW[k][l] += 0.0;
226
227 }
228 }
229
230 // tmp: verbose output
231 if (DEBUG) {
232 if (max(abs(LsusdW[k][l]),abs(RsusdW[k][l]))> 1e-6) {
233 cout << " LsusdW[" << k << "][" << l << "] = "
234 << scientific << setw(10) << LsusdW[k][l]
235 << " RsusdW[" << k << "][" << l << "] = "
236 << scientific << setw(10) << RsusdW[k][l] << endl;
237 }
238 }
239
240 }
241 }
242
243 // Now we come to the ones with really many indices
244
245 // Construct ~chi0 couplings (allow for 5 neutralinos in NMSSM)
246 for (int i=1;i<=nNeut;i++) {
247
248 // Ni1, Ni2, Ni3, Ni4, Ni5
249 complex ni1,ni2,ni3,ni4,ni5;
250 if (not isNMSSM) {
251 ni1=complex( slhaPtr->nmix(i,1), slhaPtr->imnmix(i,1) );
252 ni2=complex( slhaPtr->nmix(i,2), slhaPtr->imnmix(i,2) );
253 ni3=complex( slhaPtr->nmix(i,3), slhaPtr->imnmix(i,3) );
254 ni4=complex( slhaPtr->nmix(i,4), slhaPtr->imnmix(i,4) );
255 ni5=complex( 0.0, 0.0);
256 } else {
257 ni1=complex( slhaPtr->nmnmix(i,1), slhaPtr->imnmnmix(i,1) );
258 ni2=complex( slhaPtr->nmnmix(i,2), slhaPtr->imnmnmix(i,2) );
259 ni3=complex( slhaPtr->nmnmix(i,3), slhaPtr->imnmnmix(i,3) );
260 ni4=complex( slhaPtr->nmnmix(i,4), slhaPtr->imnmnmix(i,4) );
261 ni5=complex( slhaPtr->nmnmix(i,5), slhaPtr->imnmnmix(i,5) );
262 }
263
264 // Change to positive mass convention
265 complex iRot( 0., 1.);
266 if (slhaPtr->mass(idNeut(i)) < 0.) {
267 ni1 *= iRot;
268 ni2 *= iRot;
269 ni3 *= iRot;
270 ni4 *= iRot;
271 ni5 *= iRot;
272 }
273
274 // ~chi0 [i] ~chi0 [j] Z : loop over [j]
275 for (int j=1; j<=nNeut; j++) {
276
277 // neutralino [j] higgsino components
278 complex nj3, nj4;
279 if (not isNMSSM) {
280 nj3=complex( slhaPtr->nmix(j,3), slhaPtr->imnmix(j,3) );
281 nj4=complex( slhaPtr->nmix(j,4), slhaPtr->imnmix(j,4) );
282 } else {
283 nj3=complex( slhaPtr->nmnmix(j,3), slhaPtr->imnmnmix(j,3) );
284 nj4=complex( slhaPtr->nmnmix(j,4), slhaPtr->imnmnmix(j,4) );
285 }
286 // Change to positive mass convention
287 if (slhaPtr->mass(idNeut(j)) < 0.) {
288 nj3 *= iRot;
289 nj4 *= iRot;
290 }
291
292 // ~chi0 [i] ~chi0 [j] Z : couplings
293 OLpp[i][j] = -0.5 * ni3 * conj(nj3) + 0.5 * ni4 * conj(nj4);
294 ORpp[i][j] = 0.5 * conj(ni3) * nj3 - 0.5 * conj(ni4) * nj4;
295
296 // tmp: verbose output
297 if (DEBUG) {
298 cout << " OL'' [" << i << "][" << j << "] = "
299 << scientific << setw(10) << OLpp[i][j]
300 << " OR'' [" << i << "][" << j << "] = "
301 << scientific << setw(10) << ORpp[i][j] << endl;
302 }
303
304 }
305
306 // ~chi0 [i] ~chi+ [j] W : loop over [j]
307 for (int j=1; j<=nChar; j++) {
308
309 // Chargino mixing
310 complex uj1, uj2, vj1, vj2;
311 uj1=complex( slhaPtr->umix(j,1), slhaPtr->imumix(j,1) );
312 uj2=complex( slhaPtr->umix(j,2), slhaPtr->imumix(j,2) );
313 vj1=complex( slhaPtr->vmix(j,1), slhaPtr->imvmix(j,1) );
314 vj2=complex( slhaPtr->vmix(j,2), slhaPtr->imvmix(j,2) );
315
316 // ~chi0 [i] ~chi+ [j] W : couplings
317 OL[i][j] = -1.0/sqrt(2.0)*ni4*conj(vj2)+ni2*conj(vj1);
318 OR[i][j] = 1.0/sqrt(2.0)*conj(ni3)*uj2+conj(ni2)*uj1;
319
320 // tmp: verbose output
321 if (DEBUG) {
322 cout << " OL [" << i << "][" << j << "] = "
323 << scientific << setw(10) << OL[i][j]
324 << " OR [" << i << "][" << j << "] = "
325 << scientific << setw(10) << OR[i][j] << endl;
326 }
327 }
328
329 // Charges
330 double ed = -1.0/3.0;
331 double T3d = -0.5;
332 double eu = 2.0/3.0;
333 double T3u = 0.5;
334
335 // Loop over quark [k] generation
336 for (int k=1;k<=3;k++) {
337
338 // Set quark masses
339 // Initial guess 0,0,0,mc,mb,mt with the latter from the PDT
340 double mu = particleDataPtr->m0(2*k);
341 double md = particleDataPtr->m0(2*k-1);
342 if (k == 1) { mu=0.0 ; md=0.0; }
343 if (k == 2) { md=0.0 ; mu=0.0; }
344
345 // Compute running mass from Yukawas and vevs if possible.
346 if (slhaPtr->yd.exists() && slhaPtr->hmix.exists(3)) {
347 double ykk=slhaPtr->yd(k,k);
348 double v1=slhaPtr->hmix(3)/sqrt(1+pow(tanb,2));
349 if (ykk > 0.0) md = ykk * v1 / sqrt(2.0) ;
350 }
351 if (slhaPtr->yu.exists() && slhaPtr->hmix.exists(3)) {
352 double ykk=slhaPtr->yu(k,k);
353 double v2=slhaPtr->hmix(3)/sqrt(1.0+1.0/pow(tanb,2));
354 if (ykk > 0.0) mu = ykk * v2 / sqrt(2.0) ;
355 }
356
357 // tmp: verbose output
358 if (DEBUG) {
359 cout << " Gen = " << k << " mu = " << mu << " md = " << md
360 << " yUU,DD = " << slhaPtr->yu(k,k) << ","
361 << slhaPtr->yd(k,k) << endl;
362 }
363
364 // Loop over squark [j] flavour
365 for (int j=1;j<=6;j++) {
366
367 // Squark mixing
368 complex Rdjk = complex(Rd(j,k), imRd(j,k) );
369 complex Rdjk3 = complex(Rd(j,k+3),imRd(j,k+3));
370 complex Rujk = complex(Ru(j,k), imRu(j,k) );
371 complex Rujk3 = complex(Ru(j,k+3),imRu(j,k+3));
372
373 // ~d[j] d[k] ~chi0[i]
374 double rt2 = sqrt(2.0);
375 LsddX[j][k][i] = ((ed-T3d)*sinW/cosW*ni1 + T3d*ni2)*conj(Rdjk)/rt2
376 + md*ni3*conj(Rdjk3)/2.0/rt2/mW/cosb;
377 RsddX[j][k][i] = -ed*sinW/cosW*conj(ni1)*conj(Rdjk3)/rt2
378 + md*conj(ni3)*conj(Rdjk)/2.0/rt2/mW/cosb;
379
380 // ~u[j] u[k] ~chi0[i]
381 LsuuX[j][k][i] = ((eu-T3u)*sinW/cosW*ni1 + T3u*ni2)*conj(Rujk)/rt2
382 + mu*ni4*conj(Rujk3)/2.0/rt2/mW/sinb;
383 RsuuX[j][k][i] = -eu*sinW/cosW*conj(ni1)*conj(Rujk3)/rt2
384 + mu*conj(ni4)*conj(Rujk)/2.0/rt2/mW/sinb;
385
386 if (DEBUG) {
387 if (abs(LsddX[j][k][i]) > 1e-6) {
388 // tmp: verbose output
389 cout << " LsddX[" << j << "][" << k << "][" << i << "] = "
390 << scientific << setw(10) << LsddX[j][k][i] << endl;
391 }
392 if (abs(RsddX[j][k][i]) > 1e-6) {
393 // tmp: verbose output
394 cout << " RsddX[" << j << "][" << k << "][" << i << "] = "
395 << scientific << setw(10) << RsddX[j][k][i] << endl;
396 }
397 if (abs(LsuuX[j][k][i]) > 1e-6) {
398 // tmp: verbose output
399 cout << " LsuuX[" << j << "][" << k << "][" << i << "] = "
400 << scientific << setw(10) << LsuuX[j][k][i] << endl;
401 }
402 if (abs(RsuuX[j][k][i]) > 1e-6) {
403 // tmp: verbose output
404 cout << " RsuuX[" << j << "][" << k << "][" << i << "] = "
405 << scientific << setw(10) << RsuuX[j][k][i] << endl;
406 }
407 }
408 }
409
410 }
411
412 }
413
414 // Construct ~chi+ couplings
415 for (int i=1;i<=nChar;i++) {
416
417 // Ui1, Ui2, Vi1, Vi2
418 complex ui1,ui2,vi1,vi2;
419 ui1=complex( slhaPtr->umix(i,1), slhaPtr->imumix(i,1) );
420 ui2=complex( slhaPtr->umix(i,2), slhaPtr->imumix(i,2) );
421 vi1=complex( slhaPtr->vmix(i,1), slhaPtr->imvmix(i,1) );
422 vi2=complex( slhaPtr->vmix(i,2), slhaPtr->imvmix(i,2) );
423
424 // ~chi+ [i] ~chi- [j] Z : loop over [j]
425 for (int j=1; j<=nChar; j++) {
426
427 // Chargino mixing
428 complex uj1, uj2, vj1, vj2;
429 uj1=complex( slhaPtr->umix(j,1), slhaPtr->imumix(j,1) );
430 uj2=complex( slhaPtr->umix(j,2), slhaPtr->imumix(j,2) );
431 vj1=complex( slhaPtr->vmix(j,1), slhaPtr->imvmix(j,1) );
432 vj2=complex( slhaPtr->vmix(j,2), slhaPtr->imvmix(j,2) );
433
434 // ~chi+ [i] ~chi- [j] Z : couplings
435 OLp[i][j] = -vi1*conj(vj1) - 0.5*vi2*conj(vj2)
436 + ( (i == j) ? sin2W : 0.0);
437 ORp[i][j] = -conj(ui1)*uj1 - 0.5*conj(ui2)*uj2
438 + ( (i == j) ? sin2W : 0.0);
439
440 if (DEBUG) {
441 // tmp: verbose output
442 cout << " OL' [" << i << "][" << j << "] = "
443 << scientific << setw(10) << OLp[i][j]
444 << " OR' [" << i << "][" << j << "] = "
445 << scientific << setw(10) << ORp[i][j] << endl;
446 }
447 }
448
449 // Loop over quark [l] flavour
450 for (int l=1;l<=3;l++) {
451
452 // Set quark [l] masses
453 // Initial guess 0,0,0,mc,mb,mt with the latter from the PDT
454 double mul = particleDataPtr->m0(2*l);
455 double mdl = particleDataPtr->m0(2*l-1);
456 if (l == 1) { mul=0.0 ; mdl=0.0; }
457 if (l == 2) { mdl=0.0 ; mul=0.0; }
458
459 // Compute running mass from Yukawas and vevs if possible.
460 if (slhaPtr->yd.exists() && slhaPtr->hmix.exists(3)) {
461 double yll=slhaPtr->yd(l,l);
462 double v1=slhaPtr->hmix(3)/sqrt(1+pow(tanb,2));
463 if (yll > 0.0) mdl = yll * v1 / sqrt(2.0) ;
464 }
465 if (slhaPtr->yu.exists() && slhaPtr->hmix.exists(3)) {
466 double yll=slhaPtr->yu(l,l);
467 double v2=slhaPtr->hmix(3)/sqrt(1.0+1.0/pow(tanb,2));
468 if (yll > 0.0) mul = yll * v2 / sqrt(2.0) ;
469 }
470
471 // Loop over squark [j] flavour
472 for (int j=1;j<=6;j++) {
473
474 // Loop over off-diagonal quark [k] generation
475 for (int k=1;k<=3;k++) {
476
477 // Set quark [k] masses
478 // Initial guess 0,0,0,0,mb,mt with the latter from the PDT
479 double muk = particleDataPtr->m0(2*k);
480 double mdk = particleDataPtr->m0(2*k-1);
481 if (k == 1) { muk=0.0 ; mdk=0.0; }
482 if (k == 2) { mdk=0.0 ; muk=0.0; }
483
484 // Compute running mass from Yukawas and vevs if possible.
485 if (slhaPtr->yd.exists() && slhaPtr->hmix.exists(3)) {
486 double ykk=slhaPtr->yd(k,k);
487 double v1=slhaPtr->hmix(3)/sqrt(1+pow(tanb,2));
488 if (ykk > 0.0) mdk = ykk * v1 / sqrt(2.0) ;
489 }
490 if (slhaPtr->yu.exists() && slhaPtr->hmix.exists(3)) {
491 double ykk=slhaPtr->yu(k,k);
492 double v2=slhaPtr->hmix(3)/sqrt(1.0+1.0/pow(tanb,2));
493 if (ykk > 0.0) muk = ykk * v2 / sqrt(2.0) ;
494 }
495
496 // CKM matrix (use Pythia one if no SLHA)
497 // (NB: could also try input one if no running one found, but
498 // would then need to compute from Wolfenstein)
499 complex Vlk=coupSMPtr->VCKMgen(l,k);
500 complex Vkl=coupSMPtr->VCKMgen(k,l);
501 if (slhaPtr->vckm.exists()) {
502 Vlk=complex(slhaPtr->vckm(l,k),slhaPtr->imvckm(l,k));
503 Vkl=complex(slhaPtr->vckm(k,l),slhaPtr->imvckm(k,l));
504 }
505
506 // Squark mixing
507 complex Rdjk = complex(Rd(j,k), imRd(j,k) );
508 complex Rdjk3 = complex(Rd(j,k+3),imRd(j,k+3));
509 complex Rujk = complex(Ru(j,k), imRu(j,k) );
510 complex Rujk3 = complex(Ru(j,k+3),imRu(j,k+3));
511
512 // sqrt(2)
513 double rt2 = sqrt(2.0);
514
515 // ~d[j] u[l] ~chi+[i]
516 LsduX[j][l][i] += 0.5*(ui1*conj(Rdjk)
517 - mdk*ui2*conj(Rdjk3)/rt2/mW/cosb)*Vlk;
518 RsduX[j][l][i] -= 0.5*mul*conj(vi2)*Vlk*Rdjk/rt2/mW/sinb;
519
520 // ~u[j] d[l] ~chi+[i]
521 LsudX[j][l][i] += 0.5*(vi1*conj(Rujk)
522 - muk*vi2*conj(Rujk3)/rt2/mW/sinb)*conj(Vkl);
523 RsudX[j][l][i] -= 0.5*mdl*conj(ui2)*conj(Vkl)*conj(Rujk)/rt2/mW/cosb;
524
525 }
526
527 if (DEBUG) {
528 if (max(abs(LsduX[j][l][i]),abs(RsduX[j][l][i])) > 1e-6) {
529 // tmp: verbose output
530 cout << " LsduX[" << j << "][" << l << "][" << i << "] = "
531 << scientific << setw(10) << LsduX[j][l][i];
532 cout << " RsduX[" << j << "][" << l << "][" << i << "] = "
533 << scientific << setw(10) << RsduX[j][l][i] << endl;
534 }
535 if (max(abs(LsudX[j][l][i]),abs(RsudX[j][l][i])) > 1e-6) {
536 // tmp: verbose output
537 cout << " LsudX[" << j << "][" << l << "][" << i << "] = "
538 << scientific << setw(10) << LsudX[j][l][i];
539 cout << " RsudX[" << j << "][" << l << "][" << i << "] = "
540 << scientific << setw(10) << RsudX[j][l][i] << endl;
541 }
542 }
543
544 }
545
546
547 }
548
549 }
550
551 // SLHA2 compatibility
552 // Check whether scalar particle masses are ordered
553 bool orderedQ = true;
554 bool orderedL = true;
555 for (int j=1;j<=5;j++) {
556 if (particleDataPtr->m0(idSlep(j+1)) < particleDataPtr->m0(idSlep(j)))
557 orderedL = false;
558 if (particleDataPtr->m0(idSup(j+1)) < particleDataPtr->m0(idSup(j)))
559 orderedQ = false;
560 if (particleDataPtr->m0(idSdown(j+1)) <particleDataPtr->m0(idSdown(j)))
561 orderedQ = false;
562 }
563
564 // If ordered, change sparticle labels to mass-ordered enumeration
565 for (int i=1;i<=6;i++) {
566 ostringstream indx;
567 indx << i;
568 string uName = "~u_"+indx.str();
569 string dName = "~d_"+indx.str();
570 string lName = "~e_"+indx.str();
571 ParticleDataEntry* pdePtr;
572 if (orderedQ) {
573 pdePtr = particleDataPtr->particleDataEntryPtr(idSup(i));
574 pdePtr->setNames(uName,uName+"bar");
575 pdePtr = particleDataPtr->particleDataEntryPtr(idSdown(i));
576 pdePtr->setNames(dName,dName+"bar");
577 }
578 if (orderedL) {
579 pdePtr = particleDataPtr->particleDataEntryPtr(idSlep(i));
580 pdePtr->setNames(lName,lName+"bar");
581 }
582 }
583
584 // Let everyone know we are ready
585 isInit = true;
586
587}
588
589//--------------------------------------------------------------------------
590
591// Return neutralino flavour codes.
592
593int CoupSUSY::idNeut(int idChi) {
594 int id = 0;
595 if (idChi == 1) id = 1000022;
596 else if (idChi == 2) id = 1000023;
597 else if (idChi == 3) id = 1000025;
598 else if (idChi == 4) id = 1000035;
599 else if (idChi == 5) id = 1000045;
600 return id;
601}
602
603//--------------------------------------------------------------------------
604
605// Return chargino flavour codes.
606
607int CoupSUSY::idChar(int idChi) {
608 int id = 0;
609 if (idChi == 1) id = 1000024;
610 else if (idChi == -1) id = -1000024;
611 else if (idChi == 2) id = 1000037;
612 else if (idChi == -2) id = -1000037;
613 return id;
614}
615
616//--------------------------------------------------------------------------
617
618// Return sup flavour codes.
619
620int CoupSUSY::idSup(int iSup) {
621 int id = 0;
622 int sgn = ( iSup > 0 ) ? 1 : -1;
623 iSup = abs(iSup);
624 if (iSup == 1) id = 1000002;
625 else if (iSup == 2) id = 1000004;
626 else if (iSup == 3) id = 1000006;
627 else if (iSup == 4) id = 2000002;
628 else if (iSup == 5) id = 2000004;
629 else if (iSup == 6) id = 2000006;
630 return sgn*id;
631}
632
633//--------------------------------------------------------------------------
634
635// Return sdown flavour codes
636
637int CoupSUSY::idSdown(int iSdown) {
638 int id = 0;
639 int sgn = ( iSdown > 0 ) ? 1 : -1;
640 iSdown = abs(iSdown);
641 if (iSdown == 1) id = 1000001;
642 else if (iSdown == 2) id = 1000003;
643 else if (iSdown == 3) id = 1000005;
644 else if (iSdown == 4) id = 2000001;
645 else if (iSdown == 5) id = 2000003;
646 else if (iSdown == 6) id = 2000005;
647 return sgn*id;
648}
649
650//--------------------------------------------------------------------------
651
652// Function to return slepton flavour codes
653
654int CoupSUSY::idSlep(int iSlep) {
655 int id = 0;
656 int sgn = ( iSlep > 0 ) ? 1 : -1;
657 iSlep = abs(iSlep);
658 if (iSlep == 1) id = 1000011;
659 else if (iSlep == 2) id = 1000013;
660 else if (iSlep == 3) id = 1000015;
661 else if (iSlep == 4) id = 2000011;
662 else if (iSlep == 5) id = 2000013;
663 else if (iSlep == 6) id = 2000015;
664 return sgn*id;
665}
666
667//--------------------------------------------------------------------------
668
669// Return a particle name, given the PDG code.
670
671string CoupSUSY::getName(int pdgCode) {
672
673 // Absolute value and corresponding SM code
674 int codeA = abs(pdgCode);
675 int idSM = codeA % 1000000;
676
677 // Name
678 string name;
679
680 // Flag to indicate whether SLHA1 or SLHA2
681 bool slha1 = false;
682
683 // SM particles
684 if (codeA == idSM) {
685 // Neutrinos
686 if (!slhaPtr->upmns.exists()) slha1=true;
687 if (codeA == 12) name = (slha1) ? "nu_e" : "nu_1";
688 if (codeA == 14) name = (slha1) ? "nu_mu" : "nu_2";
689 if (codeA == 16) name = (slha1) ? "nu_tau" : "nu_3";
690 }
691
692 // Squarks
693 else if ( idSM <= 6) {
694 // up squarks
695 if (idSM % 2 == 0) {
696 // If SLHA1, return old PDG names
697 if (slhaPtr->stopmix.exists()) slha1 = true;
698 if (codeA == 1000002) name = (slha1) ? "~u_L" : "~u_1";
699 if (codeA == 1000004) name = (slha1) ? "~c_L" : "~u_2";
700 if (codeA == 1000006) name = (slha1) ? "~t_1" : "~u_3";
701 if (codeA == 2000002) name = (slha1) ? "~u_R" : "~u_4";
702 if (codeA == 2000004) name = (slha1) ? "~c_R" : "~u_5";
703 if (codeA == 2000006) name = (slha1) ? "~t_2" : "~u_6";
704 }
705 // down squarks
706 else {
707 // If SLHA1, return old PDG names
708 if (slhaPtr->sbotmix.exists()) slha1 = true;
709 if (codeA == 1000001) name = (slha1) ? "~d_L" : "~d_1";
710 if (codeA == 1000003) name = (slha1) ? "~s_L" : "~d_2";
711 if (codeA == 1000005) name = (slha1) ? "~b_1" : "~d_3";
712 if (codeA == 2000001) name = (slha1) ? "~d_R" : "~d_4";
713 if (codeA == 2000003) name = (slha1) ? "~s_R" : "~d_5";
714 if (codeA == 2000005) name = (slha1) ? "~b_2" : "~d_6";
715 }
716 if (pdgCode < 0) name += "bar";
717 }
718
719 // Sleptons
720 else if ( idSM <= 19 ) {
721
722 // Sneutrinos
723 if (idSM % 2 == 0) {
724 // If SLHA1, return old PDG names
725 if (slhaPtr->staumix.exists()) slha1 = true;
726 if (codeA == 1000012) name = (slha1) ? "~nu_eL" : "~nu_1";
727 if (codeA == 1000014) name = (slha1) ? "~nu_muL" : "~nu_2";
728 if (codeA == 1000016) name = (slha1) ? "~nu_tauL" : "~nu_3";
729 if (codeA == 2000012) name = (slha1) ? "~nu_eR" : "~nu_4";
730 if (codeA == 2000014) name = (slha1) ? "~nu_muR" : "~nu_5";
731 if (codeA == 2000016) name = (slha1) ? "~nu_tauR" : "~nu_6";
732 if (pdgCode < 0) name += "bar";
733 }
734 // charged sleptons
735 else {
736 // If SLHA1, return old PDG names
737 if (slhaPtr->staumix.exists()) slha1 = true;
738 if (codeA == 1000011) name = (slha1) ? "~e_L" : "~e_1";
739 if (codeA == 1000013) name = (slha1) ? "~mu_L" : "~e_2";
740 if (codeA == 1000015) name = (slha1) ? "~tau_1" : "~e_3";
741 if (codeA == 2000011) name = (slha1) ? "~e_R" : "~e_4";
742 if (codeA == 2000013) name = (slha1) ? "~mu_R" : "~e_5";
743 if (codeA == 2000015) name = (slha1) ? "~tau_2" : "~e_6";
744 if (pdgCode < 0) name += "-";
745 else name += "+";
746 }
747 }
748
749 else if (codeA == 1000021) name = "~g";
750 else if (codeA == 1000022) name = "~chi_10";
751 else if (codeA == 1000023) name = "~chi_20";
752 else if (codeA == 1000024) name = (pdgCode > 0) ? "~chi_1+" : "~chi_1-";
753 else if (codeA == 1000025) name = "~chi_30";
754 else if (codeA == 1000035) name = "~chi_40";
755 else if (codeA == 1000037) name = (pdgCode > 0) ? "~chi_2+" : "~chi_2-";
756
757 return name;
758
759}
760
761//==========================================================================
762
763} // end namespace Pythia8
764
765