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