StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
SusyCouplings.cc
1 // SusyCouplings.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2020 Torbjorn Sjostrand.
3 // Main authors of this file: N. Desai, P. Skands
4 // PYTHIA is licenced under the GNU GPL v2 or later, 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 // supersymmetric couplings class.
9 
10 #include "Pythia8/ParticleData.h"
11 #include "Pythia8/SusyCouplings.h"
12 
13 namespace Pythia8 {
14 
15 //==========================================================================
16 
17 // The CoupSUSY class.
18 
19 //--------------------------------------------------------------------------
20 
21 // Constants: could be changed here if desired, but normally should not.
22 // These are of technical nature, as described for each.
23 
24 //--------------------------------------------------------------------------
25 
26 // Initialize SM+SUSY couplings (only performed once).
27 
28 void CoupSUSY::initSUSY (SusyLesHouches* slhaPtrIn, Info* infoPtrIn) {
29 
30  // Save pointers.
31  infoPtr = infoPtrIn;
32  slhaPtr = slhaPtrIn;
33  settingsPtr = infoPtr->settingsPtr;
34  particleDataPtr = infoPtr->particleDataPtr;
35  coupSMPtr = infoPtr->coupSMPtr;
36 
37  // Allow verbose printout for debug purposes.
38  bool DBSUSY = settingsPtr->mode("SLHA:verbose") > 2 ? true : false;
39 
40  // Only initialize SUSY parts if SUSY is actually switched on
41  if (!slhaPtr->modsel.exists()) return;
42 
43  // Is NMSSM switched on?
44  isNMSSM = (slhaPtr->modsel(3) != 1 ? false : true);
45  settingsPtr->flag("SLHA:NMSSM",isNMSSM);
46  int nNeut = (isNMSSM ? 5 : 4);
47  int nChar = 2;
48 
49  // Initialize pole masses
50  mZpole = particleDataPtr->m0(23);
51  wZpole = particleDataPtr->mWidth(23);
52  mWpole = particleDataPtr->m0(24);
53  wWpole = particleDataPtr->mWidth(24);
54 
55  // Running masses and weak mixing angle
56  // (default to pole values if no running available)
57  mW = mWpole;
58  mZ = mZpole;
59  sin2W = coupSMPtr->sin2thetaW();
60 
61  if (settingsPtr->mode("SUSY:sin2thetaWMode") == 3) {
62  // Possibility to force on-shell sin2W definition, mostly intended for
63  // cross-checks
64  sin2W = 1.0 - pow(mW/mZ,2);
65  } else if (settingsPtr->mode("SUSY:sin2thetaWMode") == 2){
66  // Possibility to use running sin2W definition, derived from gauge
67  // couplings in running SLHA blocks (normally defined at SUSY scale)
68  if(slhaPtr->gauge.exists(1) && slhaPtr->gauge.exists(2)
69  && slhaPtr->hmix.exists(3)) {
70  double gp=slhaPtr->gauge(1);
71  double g =slhaPtr->gauge(2);
72  double v =slhaPtr->hmix(3);
73  mW = g * v / 2.0;
74  mZ = sqrt(pow(gp,2)+pow(g,2)) * v / 2.0;
75  //double tan2W = pow2(gp)/pow2(g);
76  //if (DBSUSY) cout << " tan2W = " << tan2W << endl;
77  sin2W = pow2(gp)/(pow2(g)+pow2(gp));
78  } else {
79  infoPtr->errorMsg("Warning in CoupSUSY::initSUSY: Block GAUGE"
80  " not found or incomplete; using sin(thetaW) at mZ");
81  }
82  }
83 
84  sinW = sqrt(sin2W);
85  cosW = sqrt(1.0-sin2W);
86 
87  // Tan(beta)
88  // By default, use the running one in HMIX (if not found, use MINPAR)
89 
90  if(slhaPtr->hmix.exists(2))
91  tanb = slhaPtr->hmix(2);
92  else{
93  tanb = slhaPtr->minpar(3);
94  infoPtr->errorMsg("Warning in CoupSUSY::initSUSY: Block HMIX"
95  " not found or incomplete; using MINPAR tan(beta)");
96  }
97  cosb = sqrt( 1.0 / (1.0 + tanb*tanb) );
98  sinb = sqrt( max(0.0, 1.0 - cosb*cosb));
99 
100  double beta = acos(cosb);
101 
102  // Verbose output
103  if (DBSUSY) {
104  cout << " sin2W(Q) = " << sin2W << " mW(Q) = " << mW
105  << " mZ(Q) = " << mZ << endl;
106  cout << " vev(Q) = " << slhaPtr->hmix(3) << " tanb(Q) = " << tanb
107  << endl;
108  for (int i=1;i<=3;i++) {
109  for (int j=1;j<=3;j++) {
110  cout << " VCKM [" << i << "][" << j << "] = "
111  << scientific << setw(10) << coupSMPtr->VCKMgen(i,j) << endl;
112  }
113  }
114  }
115 
116  // Higgs sector
117  if(slhaPtr->alpha.exists()) {
118  alphaHiggs = slhaPtr->alpha();
119  }
120  // If RPV, assume alpha = asin(RVHMIX(1,2)) (ignore Higgs-Sneutrino mixing)
121  else if (slhaPtr->modsel(4) == 1) {
122  alphaHiggs = asin(slhaPtr->rvhmix(1,2));
123  infoPtr->errorMsg("Info from CoupSUSY::initSUSY:","Extracting angle"
124  " alpha from RVHMIX", true);
125  } else {
126  infoPtr->errorMsg("Info from CoupSUSY::initSUSY:","Block ALPHA"
127  " not found; using alpha = beta.", true);
128  // Define approximate alpha by simple SM limit
129  alphaHiggs = atan(tanb);
130  }
131 
132  if(slhaPtr->hmix.exists(1) && slhaPtr->hmix.exists(4)){
133  muHiggs = slhaPtr->hmix(1);
134  mAHiggs = sqrt(slhaPtr->hmix(4));
135  } else if (slhaPtr->rvamix.exists()){
136  mAHiggs = particleDataPtr->m0(36);
137  muHiggs = 0.0;
138  infoPtr->errorMsg("Warning from CoupSUSY::initSUSY: Block HMIX"
139  " not found or incomplete","; setting mu = 0.");
140  } else{
141  infoPtr->errorMsg("Warning from CoupSUSY::initSUSY: Block HMIX"
142  " not found or incomplete","; setting mu = mA = 0.");
143  muHiggs = 0.0;
144  mAHiggs = 0.0;
145  }
146 
147  // Pass SLHA input to 2HDM sector
148 
149  double sba = sin(beta-alphaHiggs);
150  double cba = cos(beta-alphaHiggs);
151  double cosa = cos(alphaHiggs);
152  double sina = sin(alphaHiggs);
153 
154  // h0
155  settingsPtr->parm("HiggsH1:coup2d", -sina/cosb);
156  settingsPtr->parm("HiggsH1:coup2u", cosa/sinb);
157  settingsPtr->parm("HiggsH1:coup2l", cosa/sinb);
158  settingsPtr->parm("HiggsH1:coup2Z", sba);
159  settingsPtr->parm("HiggsH1:coup2W", sba);
160  // H0
161  settingsPtr->parm("HiggsH2:coup2d", cosa/cosb);
162  settingsPtr->parm("HiggsH2:coup2u", sina/sinb);
163  settingsPtr->parm("HiggsH2:coup2l", sina/sinb);
164  settingsPtr->parm("HiggsH2:coup2Z", cba);
165  settingsPtr->parm("HiggsH2:coup2W", cba);
166  settingsPtr->parm("HiggsH2:coup2H1Z", 0.0);
167  settingsPtr->parm("HiggsH2:coup2HchgW", sba);
168  // A0
169  settingsPtr->parm("HiggsA3:coup2d", tanb);
170  settingsPtr->parm("HiggsA3:coup2u", cosb/sinb);
171  settingsPtr->parm("HiggsA3:coup2l", cosb/sinb);
172  settingsPtr->parm("HiggsA3:coup2Z", 0.0);
173  settingsPtr->parm("HiggsA3:coup2W", 0.0);
174  settingsPtr->parm("HiggsA3:coup2H1Z", cba);
175  settingsPtr->parm("HiggsA3:coup2H2Z", sba);
176  settingsPtr->parm("HiggsA3:coup2HchgW", 1.0);
177 
178  // H^+
179  settingsPtr->parm("HiggsHchg:tanBeta", tanb);
180  settingsPtr->parm("HiggsHchg:coup2H1W", cba);
181  settingsPtr->parm("HiggsHchg:coup2H2W", sba);
182 
183  // Triple higgs couplings
184 
185  double cbpa = cos(beta+alphaHiggs);
186  double sbpa = sin(beta+alphaHiggs);
187 
188  settingsPtr->parm("HiggsH1:coup2Hchg", cos(2*beta)*sbpa + 2*pow2(cosW)*sba);
189  settingsPtr->parm("HiggsH2:coup2Hchg", -cos(2*beta)*cbpa + 2*pow2(cosW)*cba);
190  settingsPtr->parm("HiggsH2:coup2H1H1", 2*sin(2*alphaHiggs)*sbpa
191  - cos(2*alphaHiggs)*cbpa);
192  settingsPtr->parm("HiggsH2:coup2A3A3", -cos(2*beta)*cbpa);
193  settingsPtr->parm("HiggsH2:coup2A3H1", 0.0);
194  settingsPtr->parm("HiggsA3:coup2H1H1", 0.0);
195  settingsPtr->parm("HiggsA3:coup2Hchg", 0.0);
196 
197  // Shorthand for squark mixing matrices
198  LHmatrixBlock<6> Ru(slhaPtr->usqmix);
199  LHmatrixBlock<6> Rd(slhaPtr->dsqmix);
200  LHmatrixBlock<6> imRu(slhaPtr->imusqmix);
201  LHmatrixBlock<6> imRd(slhaPtr->imdsqmix);
202 
203  // Construct ~g couplings
204  for (int i=1; i<=6; i++) {
205  for (int j=1; j<=3; j++) {
206  LsddG[i][j] = complex( Rd(i,j) , imRd(i,j));
207  RsddG[i][j] = complex(-Rd(i,j+3), -imRd(i,j+3));
208  LsuuG[i][j] = complex( Ru(i,j) , imRu(i,j));
209  RsuuG[i][j] = complex(-Ru(i,j+3), -imRu(i,j+3));
210 
211  if (DBSUSY) {
212  cout << " Lsddg [" << i << "][" << j << "] = "
213  << scientific << setw(10) << LsddG[i][j]
214  << " RsddG [" << i << "][" << j << "] = "
215  << scientific << setw(10) << RsddG[i][j] << endl;
216  cout << " Lsuug [" << i << "][" << j << "] = "
217  << scientific << setw(10) << LsuuG[i][j]
218  << " RsuuG [" << i << "][" << j << "] = "
219  << scientific << setw(10) << RsuuG[i][j] << endl;
220  }
221  }
222  }
223 
224  // Construct qqZ couplings
225  for (int i=1 ; i<=6 ; i++) {
226 
227  // q[i] q[i] Z (def with extra factor 2 compared to [Okun])
228  LqqZ[i] = coupSMPtr->af(i) - 2.0*coupSMPtr->ef(i)*sin2W ;
229  RqqZ[i] = - 2.0*coupSMPtr->ef(i)*sin2W ;
230 
231  // tmp: verbose output
232  if (DBSUSY) {
233  cout << " LqqZ [" << i << "][" << i << "] = "
234  << scientific << setw(10) << LqqZ[i]
235  << " RqqZ [" << i << "][" << i << "] = "
236  << scientific << setw(10) << RqqZ[i] << endl;
237  }
238  }
239 
240  // Construct ~q~qZ couplings
241  for (int i=1 ; i<=6 ; i++) {
242 
243  // Squarks can have off-diagonal couplings as well
244  for (int j=1 ; j<=6 ; j++) {
245 
246  // ~d[i] ~d[j] Z
247  LsdsdZ[i][j] = 0.0;
248  RsdsdZ[i][j] = 0.0;
249  for (int k=1;k<=3;k++) {
250  complex Rdik = complex(Rd(i,k), imRd(i,k) );
251  complex Rdjk = complex(Rd(j,k), imRd(j,k) );
252  complex Rdik3 = complex(Rd(i,k+3),imRd(i,k+3));
253  complex Rdjk3 = complex(Rd(j,k+3),imRd(j,k+3));
254  LsdsdZ[i][j] += LqqZ[1] * (Rdik*conj(Rdjk));
255  RsdsdZ[i][j] += RqqZ[1] * (Rdik3*conj(Rdjk3));
256  }
257 
258  // ~u[i] ~u[j] Z
259  LsusuZ[i][j] = 0.0;
260  RsusuZ[i][j] = 0.0;
261  for (int k=1;k<=3;k++) {
262  complex Ruik = complex(Ru(i,k) ,imRu(i,k) );
263  complex Rujk = complex(Ru(j,k) ,imRu(j,k) );
264  complex Ruik3 = complex(Ru(i,k+3),imRu(i,k+3));
265  complex Rujk3 = complex(Ru(j,k+3),imRu(j,k+3));
266  LsusuZ[i][j] += LqqZ[2] * (Ruik*conj(Rujk));
267  RsusuZ[i][j] += RqqZ[2] * (Ruik3*conj(Rujk3));
268  }
269 
270  // tmp: verbose output
271  if (DBSUSY) {
272  if (max(norm(LsdsdZ[i][j]),norm(RsdsdZ[i][j])) > 1e-6) {
273  cout << " LsdsdZ[" << i << "][" << j << "] = "
274  << scientific << setw(10) << LsdsdZ[i][j]
275  << " RsdsdZ[" << i << "][" << j << "] = "
276  << scientific << setw(10) << RsdsdZ[i][j] << endl;
277  }
278  if (max(norm(LsusuZ[i][j]),norm(RsusuZ[i][j]))> 1e-6) {
279  cout << " LsusuZ[" << i << "][" << j << "] = "
280  << scientific << setw(10) << LsusuZ[i][j]
281  << " RsusuZ[" << i << "][" << j << "] = "
282  << scientific << setw(10) << RsusuZ[i][j] << endl;
283  }
284  }
285  }
286  }
287 
288  for(int i = 1; i < 7; i++)
289  for(int j = 1; j < 7; j++){
290  Rsl[i][j] = slhaPtr->selmix(i,j);
291  }
292 
293 
294  for(int i = 1; i < 7; i++)
295  for(int j = 1; j < 7; j++){
296  if (i < 4 && j < 4) Rsv[i][j] = slhaPtr->snumix(i,j);
297  else Rsv[i][j] = 0.0;
298  }
299 
300  // In RPV, the slepton mixing matrices include Higgs bosons
301  // Here we just extract the entries corresponding to the slepton PDG codes
302  // I.e., slepton-Higgs mixing is not supported.
303 
304  // Charged sleptons
305  if (slhaPtr->modsel(4) >= 1 && slhaPtr->rvlmix.exists()) {
306  for (int i=1; i<=6; ++i)
307  for (int j=1; j<=6; ++j)
308  Rsl[i][j] = slhaPtr->rvlmix(i+1,j+2);
309  // Check for Higgs-slepton mixing in RVLMIX
310  bool hasCrossTerms = false;
311  for (int i=2; i<=7; ++i)
312  for (int j=1; j<=2; ++j)
313  if (abs(slhaPtr->rvlmix(i,j)) > 1e-5) {
314  hasCrossTerms = true;
315  break;
316  }
317  if (hasCrossTerms)
318  infoPtr->errorMsg("Warning from CoupSUSY::initSUSY: "
319  "slepton-Higgs mixing not supported internally in PYTHIA");
320  }
321 
322  // Neutral sleptons
323  if (slhaPtr->modsel(4) >= 1 && slhaPtr->rvhmix.exists()) {
324  for (int i=1; i<=3; ++i)
325  for (int j=1; j<=3; ++j)
326  Rsv[i][j] = slhaPtr->rvhmix(i+2,j+2);
327  // Check for Higgs-sneutrino mixing in RVHMIX
328  bool hasCrossTerms = false;
329  for (int i=3; i<=7; ++i)
330  for (int j=1; j<=2; ++j)
331  if (abs(slhaPtr->rvhmix(i,j)) > 1e-5) {
332  hasCrossTerms = true;
333  break;
334  }
335  if (hasCrossTerms)
336  infoPtr->errorMsg("Warning from CoupSUSY::initSUSY: "
337  "sneutrino-Higgs mixing not supported internally in PYTHIA");
338  }
339 
340  if(DBSUSY){
341  cout<<"Rsl"<<endl;
342  for(int i=1;i<=6;i++){
343  for(int j=1;j<=6;j++){
344  cout << scientific << setw(10) << Rsl[i][j]<<" ";
345  }
346  cout<<endl;
347  }
348  cout<<"Rsv"<<endl;
349  for(int i=1;i<=6;i++){
350  for(int j=1;j<=6;j++){
351  cout << scientific << setw(10) << Rsv[i][j]<<" ";
352  }
353  cout<<endl;
354  }
355  }
356 
357 
358  // Construct llZ couplings;
359  for (int i=11 ; i<=16 ; i++) {
360 
361  LllZ[i-10] = coupSMPtr->af(i) - 2.0*coupSMPtr->ef(i)*sin2W ;
362  RllZ[i-10] = - 2.0*coupSMPtr->ef(i)*sin2W ;
363 
364  // tmp: verbose output
365  if (DBSUSY) {
366  cout << " LllZ [" << i-10 << "][" << i-10 << "] = "
367  << scientific << setw(10) << LllZ[i-10]
368  << " RllZ [" << i-10 << "][" << i-10 << "] = "
369  << scientific << setw(10) << RllZ[i-10] << endl;
370  }
371  }
372 
373  // Construct ~l~lZ couplings
374  // Initialize
375  for(int i=1;i<=6;i++){
376  for(int j=1;j<=6;j++){
377  LslslZ[i][j] = 0.0;
378  RslslZ[i][j] = 0.0;
379 
380  for (int k=1;k<=3;k++) {
381  LslslZ[i][j] += LllZ[1] * Rsl[i][k] * Rsl[j][k];
382  RslslZ[i][j] += RllZ[1] * Rsl[i][k+3] * Rsl[j][k+3];
383  }
384 
385  // ~v[i] ~v[j] Z
386  LsvsvZ[i][j] = 0.0;
387  RsvsvZ[i][j] = 0.0;
388  for (int k=1;k<=3;k++)
389  LsvsvZ[i][j] += LllZ[2] * Rsv[i][k]*Rsv[j][k];
390  }
391  }
392 
393  for(int i=1;i<=6;i++){
394  for(int j=1;j<=6;j++){
395  if (DBSUSY) {
396  if (max(norm(LsvsvZ[i][j]),norm(RsvsvZ[i][j])) > 1e-6) {
397  cout << " LsvsvZ[" << i << "][" << j << "] = "
398  << scientific << setw(10) << LsvsvZ[i][j]
399  << " RsvsvZ[" << i << "][" << j << "] = "
400  << scientific << setw(10) << RsvsvZ[i][j] << endl;
401  }
402  if (max(norm(LslslZ[i][j]),norm(RslslZ[i][j]))> 1e-6) {
403  cout << " LslslZ[" << i << "][" << j << "] = "
404  << scientific << setw(10) << LslslZ[i][j]
405  << " RslslZ[" << i << "][" << j << "] = "
406  << scientific << setw(10) << RslslZ[i][j] << endl;
407  }
408  }
409  }
410  }
411 
412  // Construct udW couplings
413  // Loop over up [i] and down [j] quark generation
414  for (int i=1;i<=3;i++) {
415  for (int j=1;j<=3;j++) {
416 
417  // CKM matrix (use Pythia one if no SLHA)
418  // (NB: could also try input one if no running one found, but
419  // would then need to compute from Wolfenstein)
420  complex Vij=coupSMPtr->VCKMgen(i,j);
421  if (slhaPtr->vckm.exists()) {
422  Vij=complex(slhaPtr->vckm(i,j),slhaPtr->imvckm(i,j));
423  }
424 
425  // u[i] d[j] W
426  LudW[i][j] = sqrt(2.0) * cosW * Vij;
427  RudW[i][j] = 0.0;
428 
429  // tmp: verbose output
430  if (DBSUSY) {
431  cout << " LudW [" << i << "][" << j << "] = "
432  << scientific << setw(10) << LudW[i][j]
433  << " RudW [" << i << "][" << j << "] = "
434  << scientific << setw(10) << RudW[i][j] << endl;
435  }
436  }
437  }
438 
439 
440  // Construct ~u~dW couplings
441  // Loop over ~u[k] and ~d[l] flavours
442  for (int k=1;k<=6;k++) {
443  for (int l=1;l<=6;l++) {
444 
445  LsusdW[k][l]=0.0;
446  RsusdW[k][l]=0.0;
447 
448  // Loop over u[i] and d[j] flavours
449  for (int i=1;i<=3;i++) {
450  for (int j=1;j<=3;j++) {
451 
452  // CKM matrix (use Pythia one if no SLHA)
453  // (NB: could also try input one if no running one found, but
454  // would then need to compute from Wolfenstein)
455  complex Vij=coupSMPtr->VCKMgen(i,j);
456  if (slhaPtr->vckm.exists()) {
457  Vij=complex(slhaPtr->vckm(i,j),slhaPtr->imvckm(i,j));
458  }
459 
460  // ~u[k] ~d[l] W (add one term for each quark flavour i,j)
461  complex Ruki = complex(Ru(k,i),imRu(k,i));
462  complex Rdlj = complex(Rd(l,j),imRd(l,j));
463  LsusdW[k][l] += sqrt(2.0) * cosW * Vij * Ruki * conj(Rdlj);
464  RsusdW[k][l] += 0.0;
465 
466  }
467  }
468 
469  // tmp: verbose output
470  if (DBSUSY) {
471  if (max(norm(LsusdW[k][l]),norm(RsusdW[k][l]))> 1e-6) {
472  cout << " LsusdW[" << k << "][" << l << "] = "
473  << scientific << setw(10) << LsusdW[k][l]
474  << " RsusdW[" << k << "][" << l << "] = "
475  << scientific << setw(10) << RsusdW[k][l] << endl;
476  }
477  }
478 
479  }
480  }
481 
482 
483 
484  // Construct lvW couplings
485  for (int i=1;i<=3;i++){
486  for (int j=1;j<=3;++j){
487  LlvW[i][j] = (i==j) ? sqrt(2.0) * cosW : 0.0 ;
488  RlvW[i][j] = 0.0;
489 
490  // tmp: verbose output
491  if (DBSUSY) {
492  cout << " LlvW [" << i << "][" << j << "] = "
493  << scientific << setw(10) << LlvW[i][j]
494  << " RlvW [" << i << "][" << j << "] = "
495  << scientific << setw(10) << RlvW[i][j] << endl;
496  }
497  }
498  }
499 
500  // Construct ~l~vW couplings
501  for (int k=1;k<=6;k++) {
502  for (int l=1;l<=6;l++) {
503  LslsvW[k][l]=0.0;
504  RslsvW[k][l]=0.0;
505 
506  if(l<=3) // Only left-handed sneutrinos
507  for(int i=1;i<=3;i++)
508  LslsvW[k][l] += sqrt(2.0) * cosW * Rsl[k][i] * Rsv[l][i];
509 
510 
511  // tmp: verbose output
512  if (DBSUSY) {
513  cout << " LslsvW [" << k << "][" << l << "] = "
514  << scientific << setw(10) << LslsvW[k][l]
515  << " RslsvW [" << k << "][" << l << "] = "
516  << scientific << setw(10) << RslsvW[k][l] << endl;
517  }
518  }
519  }
520 
521  // Now we come to the ones with really many indices
522 
523  // RPV: check for lepton-neutralino mixing (not supported)
524  if (slhaPtr->modsel(4) >= 1 && slhaPtr->rvnmix.exists()) {
525  bool hasCrossTerms = false;
526  for (int i=1; i<=3; ++i)
527  for (int j=4; j<=7; ++j)
528  if (abs(slhaPtr->rvnmix(i,j)) > 1e-5) {
529  hasCrossTerms = true;
530  break;
531  }
532  if (hasCrossTerms)
533  infoPtr->errorMsg("Warning from CoupSUSY::initSUSY: "
534  "Neutrino-Neutralino mixing not supported internally in PYTHIA");
535  }
536 
537  // Construct ~chi0 couplings (allow for 5 neutralinos in NMSSM)
538  for (int i=1;i<=nNeut;i++) {
539 
540  // Ni1, Ni2, Ni3, Ni4, Ni5
541  complex ni1,ni2,ni3,ni4,ni5;
542 
543  // In RPV, ignore neutrino-neutralino mixing
544  if (slhaPtr->modsel(4) >= 1 && slhaPtr->rvnmix.exists()) {
545  ni1=complex( slhaPtr->rvnmix(i+3,4), 0.0 );
546  ni2=complex( slhaPtr->rvnmix(i+3,5), 0.0 );
547  ni3=complex( slhaPtr->rvnmix(i+3,6), 0.0 );
548  ni4=complex( slhaPtr->rvnmix(i+3,7), 0.0 );
549  ni5=complex( 0.0, 0.0);
550  }
551  else if (!isNMSSM) {
552  ni1=complex( slhaPtr->nmix(i,1), slhaPtr->imnmix(i,1) );
553  ni2=complex( slhaPtr->nmix(i,2), slhaPtr->imnmix(i,2) );
554  ni3=complex( slhaPtr->nmix(i,3), slhaPtr->imnmix(i,3) );
555  ni4=complex( slhaPtr->nmix(i,4), slhaPtr->imnmix(i,4) );
556  ni5=complex( 0.0, 0.0);
557  } else {
558  ni1=complex( slhaPtr->nmnmix(i,1), slhaPtr->imnmnmix(i,1) );
559  ni2=complex( slhaPtr->nmnmix(i,2), slhaPtr->imnmnmix(i,2) );
560  ni3=complex( slhaPtr->nmnmix(i,3), slhaPtr->imnmnmix(i,3) );
561  ni4=complex( slhaPtr->nmnmix(i,4), slhaPtr->imnmnmix(i,4) );
562  ni5=complex( slhaPtr->nmnmix(i,5), slhaPtr->imnmnmix(i,5) );
563  }
564 
565  // Change to positive mass convention
566  complex iRot( 0., 1.);
567  if (slhaPtr->mass(idNeut(i)) < 0.) {
568  ni1 *= iRot;
569  ni2 *= iRot;
570  ni3 *= iRot;
571  ni4 *= iRot;
572  ni5 *= iRot;
573  }
574 
575  // ~chi0 [i] ~chi0 [j] Z : loop over [j]
576  for (int j=1; j<=nNeut; j++) {
577 
578  // neutralino [j] higgsino components
579  complex nj3, nj4;
580  if (slhaPtr->modsel(4) >= 1 && slhaPtr->rvnmix.exists()) {
581  nj3=complex( slhaPtr->rvnmix(i+3,6), 0.0 );
582  nj4=complex( slhaPtr->rvnmix(i+3,7), 0.0 );
583  }
584  else if (!isNMSSM) {
585  nj3=complex( slhaPtr->nmix(j,3), slhaPtr->imnmix(j,3) );
586  nj4=complex( slhaPtr->nmix(j,4), slhaPtr->imnmix(j,4) );
587  } else {
588  nj3=complex( slhaPtr->nmnmix(j,3), slhaPtr->imnmnmix(j,3) );
589  nj4=complex( slhaPtr->nmnmix(j,4), slhaPtr->imnmnmix(j,4) );
590  }
591  // Change to positive mass convention
592  if (slhaPtr->mass(idNeut(j)) < 0.) {
593  nj3 *= iRot;
594  nj4 *= iRot;
595  }
596 
597  // ~chi0 [i] ~chi0 [j] Z : couplings
598  OLpp[i][j] = -0.5 * ni3 * conj(nj3) + 0.5 * ni4 * conj(nj4);
599  ORpp[i][j] = 0.5 * conj(ni3) * nj3 - 0.5 * conj(ni4) * nj4;
600 
601  // tmp: verbose output
602  if (DBSUSY) {
603  cout << " OL'' [" << i << "][" << j << "] = "
604  << scientific << setw(10) << OLpp[i][j]
605  << " OR'' [" << i << "][" << j << "] = "
606  << scientific << setw(10) << ORpp[i][j] << endl;
607  }
608 
609  }
610 
611  // ~chi0 [i] ~chi+ [j] W : loop over [j]
612  for (int j=1; j<=nChar; j++) {
613 
614  // Chargino mixing
615  complex uj1, uj2, vj1, vj2;
616  if (slhaPtr->modsel(4)<1 || !slhaPtr->rvumix.exists()) {
617  uj1=complex( slhaPtr->umix(j,1), slhaPtr->imumix(j,1) );
618  uj2=complex( slhaPtr->umix(j,2), slhaPtr->imumix(j,2) );
619  vj1=complex( slhaPtr->vmix(j,1), slhaPtr->imvmix(j,1) );
620  vj2=complex( slhaPtr->vmix(j,2), slhaPtr->imvmix(j,2) );
621  }
622  // RPV: ignore lepton-chargino mixing
623  else {
624  uj1=complex( slhaPtr->rvumix(j+3,4), 0.0 );
625  uj2=complex( slhaPtr->rvumix(j+3,5), 0.0 );
626  vj1=complex( slhaPtr->rvvmix(j+3,4), 0.0 );
627  vj2=complex( slhaPtr->rvvmix(j+3,5), 0.0 );
628  }
629 
630  // ~chi0 [i] ~chi+ [j] W : couplings
631  OL[i][j] = -1.0/sqrt(2.0)*ni4*conj(vj2)+ni2*conj(vj1);
632  OR[i][j] = 1.0/sqrt(2.0)*conj(ni3)*uj2+conj(ni2)*uj1;
633 
634  // tmp: verbose output
635  if (DBSUSY) {
636  cout << " OL [" << i << "][" << j << "] = "
637  << scientific << setw(10) << OL[i][j]
638  << " OR [" << i << "][" << j << "] = "
639  << scientific << setw(10) << OR[i][j] << endl;
640  }
641  }
642 
643 
644  // ~qqX couplings
645  // Quark Charges
646  double ed = -1.0/3.0;
647  double T3d = -0.5;
648  double eu = 2.0/3.0;
649  double T3u = 0.5;
650 
651  // Loop over quark [k] generation
652  for (int k=1;k<=3;k++) {
653 
654  // Set quark masses
655  // Initial guess 0,0,0,mc,mb,mt with the latter from the PDT
656  double mu = particleDataPtr->m0(2*k);
657  double md = particleDataPtr->m0(2*k-1);
658  if (k == 1) { mu=0.0 ; md=0.0; }
659  if (k == 2) { md=0.0 ; mu=0.0; }
660 
661  // Compute running mass from Yukawas and vevs if possible.
662  if (slhaPtr->yd.exists() && slhaPtr->hmix.exists(3)) {
663  double ykk=slhaPtr->yd(k,k);
664  double v1=slhaPtr->hmix(3)/sqrt(1+pow(tanb,2));
665  if (ykk > 0.0) md = ykk * v1 / sqrt(2.0) ;
666  }
667  if (slhaPtr->yu.exists() && slhaPtr->hmix.exists(3)) {
668  double ykk=slhaPtr->yu(k,k);
669  double v2=slhaPtr->hmix(3)/sqrt(1.0+1.0/pow(tanb,2));
670  if (ykk > 0.0) mu = ykk * v2 / sqrt(2.0) ;
671  }
672 
673  // tmp: verbose output
674  if (DBSUSY) {
675  cout << " Gen = " << k << " mu = " << mu << " md = " << md
676  << " yUU,DD = " << slhaPtr->yu(k,k) << ","
677  << slhaPtr->yd(k,k) << endl;
678  }
679 
680  // Loop over squark [j] flavour
681  for (int j=1;j<=6;j++) {
682 
683  // Squark mixing
684  complex Rdjk = complex(Rd(j,k), imRd(j,k) );
685  complex Rdjk3 = complex(Rd(j,k+3),imRd(j,k+3));
686  complex Rujk = complex(Ru(j,k), imRu(j,k) );
687  complex Rujk3 = complex(Ru(j,k+3),imRu(j,k+3));
688 
689  // ~d[j] d[k] ~chi0[i]
690  // Changed according to new notation
691  LsddX[j][k][i] = ((ed-T3d)*sinW*ni1 + T3d*cosW*ni2)*conj(Rdjk)
692  + md*cosW*ni3*conj(Rdjk3)/2.0/mW/cosb;
693  RsddX[j][k][i] = -ed*sinW*conj(ni1)*conj(Rdjk3)
694  + md*cosW*conj(ni3)*conj(Rdjk)/2.0/mW/cosb;
695 
696  // ~u[j] u[k] ~chi0[i]
697  LsuuX[j][k][i] = ((eu-T3u)*sinW*ni1 + T3u*cosW*ni2)*conj(Rujk)
698  + mu*cosW*ni4*conj(Rujk3)/2.0/mW/sinb;
699  RsuuX[j][k][i] = -eu*sinW*conj(ni1)*conj(Rujk3)
700  + mu*cosW*conj(ni4)*conj(Rujk)/2.0/mW/sinb;
701 
702  if (DBSUSY) {
703  if (norm(LsddX[j][k][i]) > 1e-6) {
704  // tmp: verbose output
705  cout << " LsddX[" << j << "][" << k << "][" << i << "] = "
706  << scientific << setw(10) << LsddX[j][k][i] << endl;
707  }
708  if (norm(RsddX[j][k][i]) > 1e-6) {
709  // tmp: verbose output
710  cout << " RsddX[" << j << "][" << k << "][" << i << "] = "
711  << scientific << setw(10) << RsddX[j][k][i] << endl;
712  }
713  if (norm(LsuuX[j][k][i]) > 1e-6) {
714  // tmp: verbose output
715  cout << " LsuuX[" << j << "][" << k << "][" << i << "] = "
716  << scientific << setw(10) << LsuuX[j][k][i] << endl;
717  }
718  if (norm(RsuuX[j][k][i]) > 1e-6) {
719  // tmp: verbose output
720  cout << " RsuuX[" << j << "][" << k << "][" << i << "] = "
721  << scientific << setw(10) << RsuuX[j][k][i] << endl;
722  }
723  }
724  }
725  }
726 
727  // Start slepton couplings
728  // Lepton Charges
729  double el = -1.0;
730  double T3l = -0.5;
731  double ev = 0.0;
732  double T3v = 0.5;
733 
734  // Need to define lepton mass
735  for (int k=1;k<=3;k++) {
736  // Set lepton masses
737  double ml(0.0);
738  if(k==3) ml = particleDataPtr->m0(15);
739 
740  for(int j=1;j<=6;j++){
741  LsllX[j][k][i] = 0;
742  RsllX[j][k][i] = 0;
743  LsvvX[j][k][i] = 0;
744  RsvvX[j][k][i] = 0;
745 
746  // ~l[j] l[k] ~chi0[i]
747  // Changed according to new notation
748  LsllX[j][k][i] = ((el-T3l)*sinW*ni1 + T3l*cosW*ni2)*Rsl[j][k]
749  + ml*cosW*ni3/2.0/mW/cosb*Rsl[j][k+3];
750  RsllX[j][k][i] = -el*sinW*conj(ni1)*Rsl[j][k+3]
751  + ml*cosW*conj(ni3)/2.0/mW/cosb*Rsl[j][k];
752 
753  if(j<3 && j==k){
754  // No sneutrino mixing; only left handed
755  // ~v[j] v[k] ~chi0[i]
756  LsvvX[j][k][i] = ((ev-T3v)*sinW*ni1 + T3v*cosW*ni2);
757  }
758 
759  if (DBSUSY) {
760  if (norm(LsllX[j][k][i]) > 1e-6) {
761  // tmp: verbose output
762  cout << " LsllX[" << j << "][" << k << "][" << i << "] = "
763  << scientific << setw(10) << LsllX[j][k][i] << endl;
764  }
765  if (norm(RsllX[j][k][i]) > 1e-6) {
766  // tmp: verbose output
767  cout << " RsllX[" << j << "][" << k << "][" << i << "] = "
768  << scientific << setw(10) << RsllX[j][k][i] << endl;
769  }
770  if (norm(LsvvX[j][k][i]) > 1e-6) {
771  // tmp: verbose output
772  cout << " LsvvX[" << j << "][" << k << "][" << i << "] = "
773  << scientific << setw(10) << LsvvX[j][k][i] << endl;
774  }
775  }
776  }
777  }
778  }
779 
780  // RPV: check for lepton-chargino mixing (not supported)
781  if (slhaPtr->modsel(4) >= 1 && slhaPtr->rvumix.exists()) {
782  bool hasCrossTerms = false;
783  for (int i=1; i<=3; ++i)
784  for (int j=4; j<=5; ++j)
785  if (abs(slhaPtr->rvumix(i,j)) > 1e-5
786  || abs(slhaPtr->rvvmix(i,j)) > 1e-5) {
787  hasCrossTerms = true;
788  break;
789  }
790  if (hasCrossTerms)
791  infoPtr->errorMsg("Warning from CoupSUSY::initSUSY: "
792  "Lepton-Chargino mixing not supported internally in PYTHIA");
793  }
794 
795  // Construct ~chi+ couplings
796  // sqrt(2)
797  double rt2 = sqrt(2.0);
798 
799  for (int i=1;i<=nChar;i++) {
800 
801  // Ui1, Ui2, Vi1, Vi2
802  complex ui1,ui2,vi1,vi2;
803  ui1=complex( slhaPtr->umix(i,1), slhaPtr->imumix(i,1) );
804  ui2=complex( slhaPtr->umix(i,2), slhaPtr->imumix(i,2) );
805  vi1=complex( slhaPtr->vmix(i,1), slhaPtr->imvmix(i,1) );
806  vi2=complex( slhaPtr->vmix(i,2), slhaPtr->imvmix(i,2) );
807 
808  // ~chi+ [i] ~chi- [j] Z : loop over [j]
809  for (int j=1; j<=nChar; j++) {
810 
811  // Chargino mixing
812  complex uj1, uj2, vj1, vj2;
813  uj1=complex( slhaPtr->umix(j,1), slhaPtr->imumix(j,1) );
814  uj2=complex( slhaPtr->umix(j,2), slhaPtr->imumix(j,2) );
815  vj1=complex( slhaPtr->vmix(j,1), slhaPtr->imvmix(j,1) );
816  vj2=complex( slhaPtr->vmix(j,2), slhaPtr->imvmix(j,2) );
817 
818  // ~chi+ [i] ~chi- [j] Z : couplings
819  OLp[i][j] = -vi1*conj(vj1) - 0.5*vi2*conj(vj2)
820  + ( (i == j) ? sin2W : 0.0);
821  ORp[i][j] = -conj(ui1)*uj1 - 0.5*conj(ui2)*uj2
822  + ( (i == j) ? sin2W : 0.0);
823 
824  if (DBSUSY) {
825  // tmp: verbose output
826  cout << " OL' [" << i << "][" << j << "] = "
827  << scientific << setw(10) << OLp[i][j]
828  << " OR' [" << i << "][" << j << "] = "
829  << scientific << setw(10) << ORp[i][j] << endl;
830  }
831  }
832 
833  // Loop over quark [l] flavour
834  for (int l=1;l<=3;l++) {
835 
836  // Set quark [l] masses
837  // Initial guess 0,0,0,mc,mb,mt with the latter from the PDT
838  double mul = particleDataPtr->m0(2*l);
839  double mdl = particleDataPtr->m0(2*l-1);
840  if (l == 1 || l == 2) { mul=0.0 ; mdl=0.0; }
841 
842  // Compute running mass from Yukawas and vevs if possible.
843  if (slhaPtr->yd.exists() && slhaPtr->hmix.exists(3)) {
844  double yll=slhaPtr->yd(l,l);
845  double v1=slhaPtr->hmix(3)/sqrt(1+pow(tanb,2));
846  if (yll > 0.0) mdl = yll * v1 / sqrt(2.0) ;
847  }
848  if (slhaPtr->yu.exists() && slhaPtr->hmix.exists(3)) {
849  double yll=slhaPtr->yu(l,l);
850  double v2=slhaPtr->hmix(3)/sqrt(1.0+1.0/pow(tanb,2));
851  if (yll > 0.0) mul = yll * v2 / sqrt(2.0) ;
852  }
853 
854  // Loop over squark [j] flavour
855  for (int j=1;j<=6;j++) {
856 
857  //Initialise to zero
858  LsduX[j][l][i] = 0.0;
859  RsduX[j][l][i] = 0.0;
860  LsudX[j][l][i] = 0.0;
861  RsudX[j][l][i] = 0.0;
862 
863  // Loop over off-diagonal quark [k] generation
864  for (int k=1;k<=3;k++) {
865 
866  // Set quark [k] masses
867  // Initial guess 0,0,0,0,mb,mt with the latter from the PDT
868  double muk = particleDataPtr->m0(2*k);
869  double mdk = particleDataPtr->m0(2*k-1);
870  if (k == 1) { muk=0.0 ; mdk=0.0; }
871  if (k == 2) { mdk=0.0 ; muk=0.0; }
872 
873  // Compute running mass from Yukawas and vevs if possible.
874  if (slhaPtr->yd.exists() && slhaPtr->hmix.exists(3)) {
875  double ykk=slhaPtr->yd(k,k);
876  double v1=slhaPtr->hmix(3)/sqrt(1+pow(tanb,2));
877  if (ykk > 0.0) mdk = ykk * v1 / sqrt(2.0) ;
878  }
879  if (slhaPtr->yu.exists() && slhaPtr->hmix.exists(3)) {
880  double ykk=slhaPtr->yu(k,k);
881  double v2=slhaPtr->hmix(3)/sqrt(1.0+1.0/pow(tanb,2));
882  if (ykk > 0.0) muk = ykk * v2 / sqrt(2.0) ;
883  }
884 
885  // CKM matrix (use Pythia one if no SLHA)
886  // (NB: could also try input one if no running one found, but
887  // would then need to compute from Wolfenstein)
888  complex Vlk=coupSMPtr->VCKMgen(l,k);
889  complex Vkl=coupSMPtr->VCKMgen(k,l);
890  if (slhaPtr->vckm.exists()) {
891  Vlk=complex(slhaPtr->vckm(l,k),slhaPtr->imvckm(l,k));
892  Vkl=complex(slhaPtr->vckm(k,l),slhaPtr->imvckm(k,l));
893  }
894 
895  // Squark mixing
896  complex Rdjk = complex(Rd(j,k), imRd(j,k) );
897  complex Rdjk3 = complex(Rd(j,k+3),imRd(j,k+3));
898  complex Rujk = complex(Ru(j,k), imRu(j,k) );
899  complex Rujk3 = complex(Ru(j,k+3),imRu(j,k+3));
900 
901 
902  // ~d[j] u[l] ~chi+[i]
903  LsduX[j][l][i] += (ui1*conj(Rdjk)
904  - mdk*ui2*conj(Rdjk3)/rt2/mW/cosb)*Vlk;
905  RsduX[j][l][i] -= mul*conj(vi2)*Vlk*conj(Rdjk)/rt2/mW/sinb;
906 
907  // ~u[j] d[l] ~chi+[i]
908  LsudX[j][l][i] += (conj(vi1)*Rujk
909  - muk*conj(vi2)*Rujk3/rt2/mW/sinb)*Vkl;
910  RsudX[j][l][i] -= mdl*ui2*Vkl*Rujk/rt2/mW/cosb;
911 
912  }
913 
914  if (DBSUSY) {
915  if (max(norm(LsduX[j][l][i]),norm(RsduX[j][l][i])) > 1e-6) {
916  // tmp: verbose output
917  cout << " LsduX[" << j << "][" << l << "][" << i << "] = "
918  << scientific << setw(10) << LsduX[j][l][i];
919  cout << " RsduX[" << j << "][" << l << "][" << i << "] = "
920  << scientific << setw(10) << RsduX[j][l][i] << endl;
921  }
922  if (max(norm(LsudX[j][l][i]),norm(RsudX[j][l][i])) > 1e-6) {
923  // tmp: verbose output
924  cout << " LsudX[" << j << "][" << l << "][" << i << "] = "
925  << scientific << setw(10) << LsudX[j][l][i];
926  cout << " RsudX[" << j << "][" << l << "][" << i << "] = "
927  << scientific << setw(10) << RsudX[j][l][i] << endl;
928  }
929  }
930  }
931  }
932 
933  // Loop over slepton [j] flavour
934  for (int j=1;j<=6;j++) {
935  for (int k=1;k<=3;k++) {
936 
937  LslvX[j][k][i] = 0.0;
938  RslvX[j][k][i] = 0.0;
939  LsvlX[j][k][i] = 0.0;
940  RsvlX[j][k][i] = 0.0;
941 
942  // Set lepton [k] masses
943  double ml = 0.0;
944  if (k == 3) ml = particleDataPtr->m0(15);
945 
946  if(j==k || j==k+3){ // No lepton mixing
947 
948  // ~l[j] v[l] ~chi+[i]
949  LslvX[j][k][i] += ui1- ml*ui2/rt2/mW/cosb;
950  RslvX[j][k][i] -= ml*conj(vi2)/rt2/mW/sinb;
951 
952  // ~v[j] l[l] ~chi+[i]
953  if(j<=3){ // No right handed sneutrinos
954  LsvlX[j][k][i] += conj(vi1) - ml*conj(vi2)/rt2/mW/sinb;
955  }
956  }
957 
958  if (DBSUSY) {
959  if (max(norm(LslvX[j][k][i]),norm(RslvX[j][k][i])) > 1e-6) {
960  // tmp: verbose output
961  cout << " LslvX[" << j << "][" << k << "][" << i << "] = "
962  << scientific << setw(10) << LslvX[j][k][i];
963  cout << " RslvX[" << j << "][" << k << "][" << i << "] = "
964  << scientific << setw(10) << RslvX[j][k][i] << endl;
965  }
966  if (max(norm(LsvlX[j][k][i]),norm(RsvlX[j][k][i])) > 1e-6) {
967  // tmp: verbose output
968  cout << " LsvlX[" << j << "][" << k << "][" << i << "] = "
969  << scientific << setw(10) << LsvlX[j][k][i];
970  cout << " RsvlX[" << j << "][" << k << "][" << i << "] = "
971  << scientific << setw(10) << RsvlX[j][k][i] << endl;
972  }
973  }
974  }
975  }
976  }
977 
978  // Shorthand for RPV couplings
979  // The input LNV lambda couplings
980  LHtensor3Block<3> rvlle(slhaPtr->rvlamlle);
981  // The input LNV lambda' couplings
982  LHtensor3Block<3> rvlqd(slhaPtr->rvlamlqd);
983  // The input BNV lambda'' couplings
984  LHtensor3Block<3> rvudd(slhaPtr->rvlamudd);
985 
986  isLLE = false;
987  isLQD = false;
988  isUDD = false;
989 
990  // Symmetry properties
991  if(rvlle.exists()){
992  isLLE = true;
993  for(int i=1;i<=3;i++){
994  for(int j=i;j<=3;j++){
995  //lambda(i,j,k)=-lambda(j,i,k)
996  for(int k=1;k<=3;k++){
997  if(i==j){
998  rvLLE[i][j][k] = 0.0;
999  }else{
1000  rvLLE[i][j][k] = rvlle(i,j,k);
1001  rvLLE[j][i][k] = -rvlle(i,j,k);
1002  }
1003  }
1004  }
1005  }
1006  }
1007 
1008  if(rvlqd.exists()){
1009  isLQD=true;
1010  for(int i=1;i<=3;i++){
1011  for(int j=1;j<=3;j++){
1012  for(int k=1;k<=3;k++){
1013  rvLQD[i][j][k] = rvlqd(i,j,k);
1014  }
1015  }
1016  }
1017  }
1018 
1019  //lambda''(k,j,i)=-lambda''(k,i,j)
1020 
1021  if(rvudd.exists()){
1022  isUDD = true;
1023  for(int i=1;i<=3;i++){
1024  for(int j=i;j<=3;j++){
1025  for(int k=1;k<=3;k++){
1026  if(i==j){
1027  rvUDD[k][i][j] = 0.0;
1028  }else{
1029  rvUDD[k][i][j] = rvudd(k,i,j);
1030  rvUDD[k][j][i] = -rvudd(k,i,j);
1031  }
1032  }
1033  }
1034  }
1035  }
1036 
1037  if(DBSUSY){
1038  for(int i=1;i<=3;i++){
1039  for(int j=1;j<=3;j++){
1040  for(int k=1;k<=3;k++){
1041  if(rvlle.exists())
1042  cout<<"LambdaLLE["<<i<<"]["<<j<<"]["<<k<<"]="<<rvLLE[i][j][k]<<" ";
1043  if(rvlqd.exists())
1044  cout<<"LambdaLQD["<<i<<"]["<<j<<"]["<<k<<"]="<<rvLQD[i][j][k]<<" ";
1045  if(rvudd.exists())
1046  cout<<"LambdaUDD["<<i<<"]["<<j<<"]["<<k<<"]="<<rvUDD[i][j][k]
1047  <<"\n";
1048  }
1049  }
1050  }
1051  }
1052 
1053  // Store the squark mixing matrix
1054  for(int i=1;i<=6;i++){
1055  for(int j=1;j<=3;j++){
1056  Rusq[i][j] = complex(Ru(i,j), imRu(i,j) );
1057  Rusq[i][j+3] = complex(Ru(i,j+3),imRu(i,j+3));
1058  Rdsq[i][j] = complex(Rd(i,j), imRd(i,j) );
1059  Rdsq[i][j+3] = complex(Rd(i,j+3),imRd(i,j+3));
1060  }
1061  }
1062 
1063  if(DBSUSY){
1064  cout<<"Ru"<<endl;
1065  for(int i=1;i<=6;i++){
1066  for(int j=1;j<=6;j++){
1067  cout << real(Rusq[i][j])<<" ";
1068  }
1069  cout<<endl;
1070  }
1071  cout<<"Rd"<<endl;
1072  for(int i=1;i<=6;i++){
1073  for(int j=1;j<=6;j++){
1074  cout << real(Rdsq[i][j])<<" ";
1075  }
1076  cout<<endl;
1077  }
1078  }
1079 
1080  // Let everyone know we are ready
1081  isInit = true;
1082 }
1083 
1084 //--------------------------------------------------------------------------
1085 
1086 // Return neutralino flavour codes.
1087 
1088 int CoupSUSY::idNeut(int idChi) {
1089 
1090  int id = 0;
1091  if (idChi == 1) id = 1000022;
1092  else if (idChi == 2) id = 1000023;
1093  else if (idChi == 3) id = 1000025;
1094  else if (idChi == 4) id = 1000035;
1095  else if (idChi == 5) id = 1000045;
1096  return id;
1097 }
1098 
1099 //--------------------------------------------------------------------------
1100 
1101 // Return chargino flavour codes.
1102 
1103 int CoupSUSY::idChar(int idChi) {
1104 
1105  int id = 0;
1106  if (idChi == 1) id = 1000024;
1107  else if (idChi == -1) id = -1000024;
1108  else if (idChi == 2) id = 1000037;
1109  else if (idChi == -2) id = -1000037;
1110  return id;
1111 }
1112 
1113 //--------------------------------------------------------------------------
1114 
1115 // Return sup flavour codes.
1116 
1117 int CoupSUSY::idSup(int iSup) {
1118 
1119  int id = 0;
1120  int sgn = ( iSup > 0 ) ? 1 : -1;
1121  iSup = abs(iSup);
1122  if (iSup == 1) id = 1000002;
1123  else if (iSup == 2) id = 1000004;
1124  else if (iSup == 3) id = 1000006;
1125  else if (iSup == 4) id = 2000002;
1126  else if (iSup == 5) id = 2000004;
1127  else if (iSup == 6) id = 2000006;
1128  return sgn*id;
1129 }
1130 
1131 //--------------------------------------------------------------------------
1132 
1133 // Return sdown flavour codes
1134 
1135 int CoupSUSY::idSdown(int iSdown) {
1136 
1137  int id = 0;
1138  int sgn = ( iSdown > 0 ) ? 1 : -1;
1139  iSdown = abs(iSdown);
1140  if (iSdown == 1) id = 1000001;
1141  else if (iSdown == 2) id = 1000003;
1142  else if (iSdown == 3) id = 1000005;
1143  else if (iSdown == 4) id = 2000001;
1144  else if (iSdown == 5) id = 2000003;
1145  else if (iSdown == 6) id = 2000005;
1146  return sgn*id;
1147 }
1148 
1149 //--------------------------------------------------------------------------
1150 
1151 // Function to return slepton flavour codes
1152 
1153 int CoupSUSY::idSlep(int iSlep) {
1154 
1155  int id = 0;
1156  int sgn = ( iSlep > 0 ) ? 1 : -1;
1157  iSlep = abs(iSlep);
1158  if (iSlep == 1) id = 1000011;
1159  else if (iSlep == 2) id = 1000013;
1160  else if (iSlep == 3) id = 1000015;
1161  else if (iSlep == 4) id = 2000011;
1162  else if (iSlep == 5) id = 2000013;
1163  else if (iSlep == 6) id = 2000015;
1164  return sgn*id;
1165 }
1166 
1167 //--------------------------------------------------------------------------
1168 
1169 // Return neutralino code; zero if not a (recognized) neutralino
1170 
1171 int CoupSUSY::typeNeut(int idPDG) {
1172  int type = 0;
1173  int idAbs = abs(idPDG);
1174  if(idAbs == 1000022) type = 1;
1175  else if(idAbs == 1000023) type = 2;
1176  else if(idAbs == 1000025) type = 3;
1177  else if(idAbs == 1000035) type = 4;
1178  else if(isNMSSM && idAbs == 1000045) type = 5;
1179  return type;
1180 
1181 }
1182 
1183 //--------------------------------------------------------------------------
1184 
1185 // Check whether particle is a Chargino
1186 
1187 int CoupSUSY::typeChar(int idPDG) {
1188  int type = 0;
1189  if(abs(idPDG) == 1000024) type = 1;
1190  else if (abs(idPDG) == 1000037)type = 2;
1191  return type;
1192 }
1193 
1194 //==========================================================================
1195 
1196 } // end namespace Pythia8
Generated by   doxygen 1.8.5