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) 2014 Torbjorn Sjostrand.
3 // Main authors of this file: N. Desai, P. Skands
4 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
5 // Please respect the MCnet Guidelines, see GUIDELINES for details.
6 
7 // Function definitions (not found in the header) for the
8 // 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  LHmatrixBlock<6> Rsl(slhaPtr->selmix);
292  LHmatrixBlock<3> Rsv(slhaPtr->snumix);
293 
294  // In RPV, the slepton mixing matrices include Higgs bosons
295  // Here we just extract the entries corresponding to the slepton PDG codes
296  // I.e., slepton-Higgs mixing is not supported.
297 
298  // Charged sleptons
299  if (slhaPtr->modsel(4) >= 1 && slhaPtr->rvlmix.exists()) {
300  for (int i=1; i<=6; ++i)
301  for (int j=1; j<=6; ++j)
302  Rsl.set(i,j,slhaPtr->rvlmix(i+1,j+2));
303  // Check for Higgs-slepton mixing in RVLMIX
304  bool hasCrossTerms = false;
305  for (int i=2; i<=7; ++i)
306  for (int j=1; j<=2; ++j)
307  if (abs(slhaPtr->rvlmix(i,j)) > 1e-5) {
308  hasCrossTerms = true;
309  break;
310  }
311  if (hasCrossTerms)
312  infoPtr->errorMsg("Warning from CoupSUSY::initSUSY: "
313  "slepton-Higgs mixing not supported in PYTHIA");
314  }
315 
316  // Neutral sleptons
317  if (slhaPtr->modsel(4) >= 1 && slhaPtr->rvhmix.exists()) {
318  for (int i=1; i<=3; ++i)
319  for (int j=1; j<=3; ++j)
320  Rsv.set(i,j,slhaPtr->rvhmix(i+2,j+2));
321  // Check for Higgs-sneutrino mixing in RVHMIX
322  bool hasCrossTerms = false;
323  for (int i=3; i<=7; ++i)
324  for (int j=1; j<=2; ++j)
325  if (abs(slhaPtr->rvhmix(i,j)) > 1e-5) {
326  hasCrossTerms = true;
327  break;
328  }
329  if (hasCrossTerms)
330  infoPtr->errorMsg("Warning from CoupSUSY::initSUSY: "
331  "sneutrino-Higgs mixing not supported in PYTHIA");
332  }
333 
334  if(DBSUSY){
335  cout<<"Rsl"<<endl;
336  for(int i=1;i<=6;i++){
337  for(int j=1;j<=6;j++){
338  cout << scientific << setw(10) << Rsl(i,j)<<" ";
339  }
340  cout<<endl;
341  }
342  cout<<"Rsv"<<endl;
343  for(int i=1;i<=6;i++){
344  for(int j=1;j<=6;j++){
345  cout << scientific << setw(10) << Rsv(i,j)<<" ";
346  }
347  cout<<endl;
348  }
349  }
350 
351 
352  // Construct llZ couplings;
353  for (int i=11 ; i<=16 ; i++) {
354 
355  LllZ[i-10] = af(i) - 2.0*ef(i)*sin2W ;
356  RllZ[i-10] = - 2.0*ef(i)*sin2W ;
357 
358  // tmp: verbose output
359  if (DBSUSY) {
360  cout << " LllZ [" << i-10 << "][" << i-10 << "] = "
361  << scientific << setw(10) << LllZ[i-10]
362  << " RllZ [" << i-10 << "][" << i-10 << "] = "
363  << scientific << setw(10) << RllZ[i-10] << endl;
364  }
365  }
366 
367  // Construct ~l~lZ couplings
368  // Initialize
369  for(int i=1;i<=6;i++){
370  for(int j=1;j<=6;j++){
371  LslslZ[i][j] = 0.0;
372  RslslZ[i][j] = 0.0;
373 
374  for (int k=1;k<=3;k++) {
375  LslslZ[i][j] += LllZ[1] * Rsl(i,k) * Rsl(j,k);
376  RslslZ[i][j] += RllZ[1] * Rsl(i,k+3) * Rsl(j,k+3);
377  }
378 
379  // ~v[i] ~v[j] Z
380  LsvsvZ[i][j] = 0.0;
381  RsvsvZ[i][j] = 0.0;
382  for (int k=1;k<=3;k++)
383  LsvsvZ[i][j] += LllZ[2] * Rsv(i,k)*Rsv(j,k);
384  }
385  }
386 
387  for(int i=1;i<=6;i++){
388  for(int j=1;j<=6;j++){
389  if (DBSUSY) {
390  if (max(abs(LsvsvZ[i][j]),abs(RsvsvZ[i][j])) > 1e-6) {
391  cout << " LsvsvZ[" << i << "][" << j << "] = "
392  << scientific << setw(10) << LsvsvZ[i][j]
393  << " RsvsvZ[" << i << "][" << j << "] = "
394  << scientific << setw(10) << RsvsvZ[i][j] << endl;
395  }
396  if (max(abs(LslslZ[i][j]),abs(RslslZ[i][j]))> 1e-6) {
397  cout << " LslslZ[" << i << "][" << j << "] = "
398  << scientific << setw(10) << LslslZ[i][j]
399  << " RslslZ[" << i << "][" << j << "] = "
400  << scientific << setw(10) << RslslZ[i][j] << endl;
401  }
402  }
403  }
404  }
405 
406  // Construct udW couplings
407  // Loop over up [i] and down [j] quark generation
408  for (int i=1;i<=3;i++) {
409  for (int j=1;j<=3;j++) {
410 
411  // CKM matrix (use Pythia one if no SLHA)
412  // (NB: could also try input one if no running one found, but
413  // would then need to compute from Wolfenstein)
414  complex Vij=VCKMgen(i,j);
415  if (slhaPtr->vckm.exists()) {
416  Vij=complex(slhaPtr->vckm(i,j),slhaPtr->imvckm(i,j));
417  }
418 
419  // u[i] d[j] W
420  LudW[i][j] = sqrt(2.0) * cosW * Vij;
421  RudW[i][j] = 0.0;
422 
423  // tmp: verbose output
424  if (DBSUSY) {
425  cout << " LudW [" << i << "][" << j << "] = "
426  << scientific << setw(10) << LudW[i][j]
427  << " RudW [" << i << "][" << j << "] = "
428  << scientific << setw(10) << RudW[i][j] << endl;
429  }
430  }
431  }
432 
433 
434  // Construct ~u~dW couplings
435  // Loop over ~u[k] and ~d[l] flavours
436  for (int k=1;k<=6;k++) {
437  for (int l=1;l<=6;l++) {
438 
439  LsusdW[k][l]=0.0;
440  RsusdW[k][l]=0.0;
441 
442  // Loop over u[i] and d[j] flavours
443  for (int i=1;i<=3;i++) {
444  for (int j=1;j<=3;j++) {
445 
446  // CKM matrix (use Pythia one if no SLHA)
447  // (NB: could also try input one if no running one found, but
448  // would then need to compute from Wolfenstein)
449  complex Vij=VCKMgen(i,j);
450  if (slhaPtr->vckm.exists()) {
451  Vij=complex(slhaPtr->vckm(i,j),slhaPtr->imvckm(i,j));
452  }
453 
454  // ~u[k] ~d[l] W (add one term for each quark flavour i,j)
455  complex Ruki = complex(Ru(k,i),imRu(k,i));
456  complex Rdlj = complex(Rd(l,j),imRd(l,j));
457  LsusdW[k][l] += sqrt(2.0) * cosW * Vij * Ruki * conj(Rdlj);
458  RsusdW[k][l] += 0.0;
459 
460  }
461  }
462 
463  // tmp: verbose output
464  if (DBSUSY) {
465  if (max(abs(LsusdW[k][l]),abs(RsusdW[k][l]))> 1e-6) {
466  cout << " LsusdW[" << k << "][" << l << "] = "
467  << scientific << setw(10) << LsusdW[k][l]
468  << " RsusdW[" << k << "][" << l << "] = "
469  << scientific << setw(10) << RsusdW[k][l] << endl;
470  }
471  }
472 
473  }
474  }
475 
476 
477 
478  // Construct lvW couplings
479  for (int i=1;i<=3;i++){
480  for (int j=1;j<=3;++j){
481  LlvW[i][j] = (i==j) ? sqrt(2.0) * cosW : 0.0 ;
482  RlvW[i][j] = 0.0;
483 
484  // tmp: verbose output
485  if (DBSUSY) {
486  cout << " LlvW [" << i << "][" << j << "] = "
487  << scientific << setw(10) << LlvW[i][j]
488  << " RlvW [" << i << "][" << j << "] = "
489  << scientific << setw(10) << RlvW[i][j] << endl;
490  }
491  }
492  }
493 
494  // Construct ~l~vW couplings
495  for (int k=1;k<=6;k++) {
496  for (int l=1;l<=6;l++) {
497  LslsvW[k][l]=0.0;
498  RslsvW[k][l]=0.0;
499 
500  if(l<=3) // Only left-handed sneutrinos
501  for(int i=1;i<=3;i++)
502  LslsvW[k][l] += sqrt(2.0) * cosW * Rsl(k,i) * Rsv(l,i);
503 
504 
505  // tmp: verbose output
506  if (DBSUSY) {
507  cout << " LslsvW [" << k << "][" << l << "] = "
508  << scientific << setw(10) << LslsvW[k][l]
509  << " RslsvW [" << k << "][" << l << "] = "
510  << scientific << setw(10) << RslsvW[k][l] << endl;
511  }
512  }
513  }
514 
515  // Now we come to the ones with really many indices
516 
517  // RPV: check for lepton-neutralino mixing (not supported)
518  if (slhaPtr->modsel(4) >= 1 && slhaPtr->rvnmix.exists()) {
519  bool hasCrossTerms = false;
520  for (int i=1; i<=3; ++i)
521  for (int j=4; j<=7; ++j)
522  if (abs(slhaPtr->rvnmix(i,j)) > 1e-5) {
523  hasCrossTerms = true;
524  break;
525  }
526  if (hasCrossTerms)
527  infoPtr->errorMsg("Warning from CoupSUSY::initSUSY: "
528  "Neutrino-Neutralino mixing not supported in PYTHIA");
529  }
530 
531  // Construct ~chi0 couplings (allow for 5 neutralinos in NMSSM)
532  for (int i=1;i<=nNeut;i++) {
533 
534  // Ni1, Ni2, Ni3, Ni4, Ni5
535  complex ni1,ni2,ni3,ni4,ni5;
536 
537  // In RPV, ignore neutrino-neutralino mixing
538  if (slhaPtr->modsel(4) >= 1 && slhaPtr->rvnmix.exists()) {
539  ni1=complex( slhaPtr->rvnmix(i+3,4), 0.0 );
540  ni2=complex( slhaPtr->rvnmix(i+3,5), 0.0 );
541  ni3=complex( slhaPtr->rvnmix(i+3,6), 0.0 );
542  ni4=complex( slhaPtr->rvnmix(i+3,7), 0.0 );
543  ni5=complex( 0.0, 0.0);
544  }
545  else if (!isNMSSM) {
546  ni1=complex( slhaPtr->nmix(i,1), slhaPtr->imnmix(i,1) );
547  ni2=complex( slhaPtr->nmix(i,2), slhaPtr->imnmix(i,2) );
548  ni3=complex( slhaPtr->nmix(i,3), slhaPtr->imnmix(i,3) );
549  ni4=complex( slhaPtr->nmix(i,4), slhaPtr->imnmix(i,4) );
550  ni5=complex( 0.0, 0.0);
551  } else {
552  ni1=complex( slhaPtr->nmnmix(i,1), slhaPtr->imnmnmix(i,1) );
553  ni2=complex( slhaPtr->nmnmix(i,2), slhaPtr->imnmnmix(i,2) );
554  ni3=complex( slhaPtr->nmnmix(i,3), slhaPtr->imnmnmix(i,3) );
555  ni4=complex( slhaPtr->nmnmix(i,4), slhaPtr->imnmnmix(i,4) );
556  ni5=complex( slhaPtr->nmnmix(i,5), slhaPtr->imnmnmix(i,5) );
557  }
558 
559  // Change to positive mass convention
560  complex iRot( 0., 1.);
561  if (slhaPtr->mass(idNeut(i)) < 0.) {
562  ni1 *= iRot;
563  ni2 *= iRot;
564  ni3 *= iRot;
565  ni4 *= iRot;
566  ni5 *= iRot;
567  }
568 
569  // ~chi0 [i] ~chi0 [j] Z : loop over [j]
570  for (int j=1; j<=nNeut; j++) {
571 
572  // neutralino [j] higgsino components
573  complex nj3, nj4;
574  if (slhaPtr->modsel(4) >= 1 && slhaPtr->rvnmix.exists()) {
575  nj3=complex( slhaPtr->rvnmix(i+3,6), 0.0 );
576  nj4=complex( slhaPtr->rvnmix(i+3,7), 0.0 );
577  }
578  else if (!isNMSSM) {
579  nj3=complex( slhaPtr->nmix(j,3), slhaPtr->imnmix(j,3) );
580  nj4=complex( slhaPtr->nmix(j,4), slhaPtr->imnmix(j,4) );
581  } else {
582  nj3=complex( slhaPtr->nmnmix(j,3), slhaPtr->imnmnmix(j,3) );
583  nj4=complex( slhaPtr->nmnmix(j,4), slhaPtr->imnmnmix(j,4) );
584  }
585  // Change to positive mass convention
586  if (slhaPtr->mass(idNeut(j)) < 0.) {
587  nj3 *= iRot;
588  nj4 *= iRot;
589  }
590 
591  // ~chi0 [i] ~chi0 [j] Z : couplings
592  OLpp[i][j] = -0.5 * ni3 * conj(nj3) + 0.5 * ni4 * conj(nj4);
593  ORpp[i][j] = 0.5 * conj(ni3) * nj3 - 0.5 * conj(ni4) * nj4;
594 
595  // tmp: verbose output
596  if (DBSUSY) {
597  cout << " OL'' [" << i << "][" << j << "] = "
598  << scientific << setw(10) << OLpp[i][j]
599  << " OR'' [" << i << "][" << j << "] = "
600  << scientific << setw(10) << ORpp[i][j] << endl;
601  }
602 
603  }
604 
605  // ~chi0 [i] ~chi+ [j] W : loop over [j]
606  for (int j=1; j<=nChar; j++) {
607 
608  // Chargino mixing
609  complex uj1, uj2, vj1, vj2;
610  if (slhaPtr->modsel(4)<1 || !slhaPtr->rvumix.exists()) {
611  uj1=complex( slhaPtr->umix(j,1), slhaPtr->imumix(j,1) );
612  uj2=complex( slhaPtr->umix(j,2), slhaPtr->imumix(j,2) );
613  vj1=complex( slhaPtr->vmix(j,1), slhaPtr->imvmix(j,1) );
614  vj2=complex( slhaPtr->vmix(j,2), slhaPtr->imvmix(j,2) );
615  }
616  // RPV: ignore lepton-chargino mixing
617  else {
618  uj1=complex( slhaPtr->rvumix(j+3,4), 0.0 );
619  uj2=complex( slhaPtr->rvumix(j+3,5), 0.0 );
620  vj1=complex( slhaPtr->rvvmix(j+3,4), 0.0 );
621  vj2=complex( slhaPtr->rvvmix(j+3,5), 0.0 );
622  }
623 
624  // ~chi0 [i] ~chi+ [j] W : couplings
625  OL[i][j] = -1.0/sqrt(2.0)*ni4*conj(vj2)+ni2*conj(vj1);
626  OR[i][j] = 1.0/sqrt(2.0)*conj(ni3)*uj2+conj(ni2)*uj1;
627 
628  // tmp: verbose output
629  if (DBSUSY) {
630  cout << " OL [" << i << "][" << j << "] = "
631  << scientific << setw(10) << OL[i][j]
632  << " OR [" << i << "][" << j << "] = "
633  << scientific << setw(10) << OR[i][j] << endl;
634  }
635  }
636 
637 
638  // ~qqX couplings
639  // Quark Charges
640  double ed = -1.0/3.0;
641  double T3d = -0.5;
642  double eu = 2.0/3.0;
643  double T3u = 0.5;
644 
645  // Loop over quark [k] generation
646  for (int k=1;k<=3;k++) {
647 
648  // Set quark masses
649  // Initial guess 0,0,0,mc,mb,mt with the latter from the PDT
650  double mu = particleDataPtr->m0(2*k);
651  double md = particleDataPtr->m0(2*k-1);
652  if (k == 1) { mu=0.0 ; md=0.0; }
653  if (k == 2) { md=0.0 ; mu=0.0; }
654 
655  // Compute running mass from Yukawas and vevs if possible.
656  if (slhaPtr->yd.exists() && slhaPtr->hmix.exists(3)) {
657  double ykk=slhaPtr->yd(k,k);
658  double v1=slhaPtr->hmix(3)/sqrt(1+pow(tanb,2));
659  if (ykk > 0.0) md = ykk * v1 / sqrt(2.0) ;
660  }
661  if (slhaPtr->yu.exists() && slhaPtr->hmix.exists(3)) {
662  double ykk=slhaPtr->yu(k,k);
663  double v2=slhaPtr->hmix(3)/sqrt(1.0+1.0/pow(tanb,2));
664  if (ykk > 0.0) mu = ykk * v2 / sqrt(2.0) ;
665  }
666 
667  // tmp: verbose output
668  if (DBSUSY) {
669  cout << " Gen = " << k << " mu = " << mu << " md = " << md
670  << " yUU,DD = " << slhaPtr->yu(k,k) << ","
671  << slhaPtr->yd(k,k) << endl;
672  }
673 
674  // Loop over squark [j] flavour
675  for (int j=1;j<=6;j++) {
676 
677  // Squark mixing
678  complex Rdjk = complex(Rd(j,k), imRd(j,k) );
679  complex Rdjk3 = complex(Rd(j,k+3),imRd(j,k+3));
680  complex Rujk = complex(Ru(j,k), imRu(j,k) );
681  complex Rujk3 = complex(Ru(j,k+3),imRu(j,k+3));
682 
683  // ~d[j] d[k] ~chi0[i]
684  // Changed according to new notation
685  LsddX[j][k][i] = ((ed-T3d)*sinW*ni1 + T3d*cosW*ni2)*conj(Rdjk)
686  + md*cosW*ni3*conj(Rdjk3)/2.0/mW/cosb;
687  RsddX[j][k][i] = -ed*sinW*conj(ni1)*conj(Rdjk3)
688  + md*cosW*conj(ni3)*conj(Rdjk)/2.0/mW/cosb;
689 
690  // ~u[j] u[k] ~chi0[i]
691  LsuuX[j][k][i] = ((eu-T3u)*sinW*ni1 + T3u*cosW*ni2)*conj(Rujk)
692  + mu*cosW*ni4*conj(Rujk3)/2.0/mW/sinb;
693  RsuuX[j][k][i] = -eu*sinW*conj(ni1)*conj(Rujk3)
694  + mu*cosW*conj(ni4)*conj(Rujk)/2.0/mW/sinb;
695 
696  if (DBSUSY) {
697  if (abs(LsddX[j][k][i]) > 1e-6) {
698  // tmp: verbose output
699  cout << " LsddX[" << j << "][" << k << "][" << i << "] = "
700  << scientific << setw(10) << LsddX[j][k][i] << endl;
701  }
702  if (abs(RsddX[j][k][i]) > 1e-6) {
703  // tmp: verbose output
704  cout << " RsddX[" << j << "][" << k << "][" << i << "] = "
705  << scientific << setw(10) << RsddX[j][k][i] << endl;
706  }
707  if (abs(LsuuX[j][k][i]) > 1e-6) {
708  // tmp: verbose output
709  cout << " LsuuX[" << j << "][" << k << "][" << i << "] = "
710  << scientific << setw(10) << LsuuX[j][k][i] << endl;
711  }
712  if (abs(RsuuX[j][k][i]) > 1e-6) {
713  // tmp: verbose output
714  cout << " RsuuX[" << j << "][" << k << "][" << i << "] = "
715  << scientific << setw(10) << RsuuX[j][k][i] << endl;
716  }
717  }
718  }
719  }
720 
721  // Start slepton couplings
722  // Lepton Charges
723  double el = -1.0;
724  double T3l = -0.5;
725  double ev = 0.0;
726  double T3v = 0.5;
727 
728  // Need to define lepton mass
729  for (int k=1;k<=3;k++) {
730  // Set lepton masses
731  double ml(0.0);
732  if(k==3) ml = particleDataPtr->m0(15);
733 
734  for(int j=1;j<=6;j++){
735  LsllX[j][k][i] = 0;
736  RsllX[j][k][i] = 0;
737  LsvvX[j][k][i] = 0;
738  RsvvX[j][k][i] = 0;
739 
740  // ~l[j] l[k] ~chi0[i]
741  // Changed according to new notation
742  LsllX[j][k][i] = ((el-T3l)*sinW*ni1 + T3l*cosW*ni2)*Rsl(j,k)
743  + ml*cosW*ni3/2.0/mW/cosb*Rsl(j,k+3);
744  RsllX[j][k][i] = -el*sinW*conj(ni1)*Rsl(j,k+3)
745  + ml*cosW*conj(ni3)/2.0/mW/cosb*Rsl(j,k);
746 
747  if(j<3 && j==k){
748  // No sneutrino mixing; only left handed
749  // ~v[j] v[k] ~chi0[i]
750  LsvvX[j][k][i] = ((ev-T3v)*sinW*ni1 + T3v*cosW*ni2);
751  }
752 
753  if (DBSUSY) {
754  if (abs(LsllX[j][k][i]) > 1e-6) {
755  // tmp: verbose output
756  cout << " LsllX[" << j << "][" << k << "][" << i << "] = "
757  << scientific << setw(10) << LsllX[j][k][i] << endl;
758  }
759  if (abs(RsllX[j][k][i]) > 1e-6) {
760  // tmp: verbose output
761  cout << " RsllX[" << j << "][" << k << "][" << i << "] = "
762  << scientific << setw(10) << RsllX[j][k][i] << endl;
763  }
764  if (abs(LsvvX[j][k][i]) > 1e-6) {
765  // tmp: verbose output
766  cout << " LsvvX[" << j << "][" << k << "][" << i << "] = "
767  << scientific << setw(10) << LsvvX[j][k][i] << endl;
768  }
769  }
770  }
771  }
772  }
773 
774  // RPV: check for lepton-chargino mixing (not supported)
775  if (slhaPtr->modsel(4) >= 1 && slhaPtr->rvumix.exists()) {
776  bool hasCrossTerms = false;
777  for (int i=1; i<=3; ++i)
778  for (int j=4; j<=5; ++j)
779  if (abs(slhaPtr->rvumix(i,j)) > 1e-5
780  || abs(slhaPtr->rvvmix(i,j)) > 1e-5) {
781  hasCrossTerms = true;
782  break;
783  }
784  if (hasCrossTerms)
785  infoPtr->errorMsg("Warning from CoupSUSY::initSUSY: "
786  "Lepton-Chargino mixing not supported in PYTHIA");
787  }
788 
789  // Construct ~chi+ couplings
790  // sqrt(2)
791  double rt2 = sqrt(2.0);
792 
793  for (int i=1;i<=nChar;i++) {
794 
795  // Ui1, Ui2, Vi1, Vi2
796  complex ui1,ui2,vi1,vi2;
797  ui1=complex( slhaPtr->umix(i,1), slhaPtr->imumix(i,1) );
798  ui2=complex( slhaPtr->umix(i,2), slhaPtr->imumix(i,2) );
799  vi1=complex( slhaPtr->vmix(i,1), slhaPtr->imvmix(i,1) );
800  vi2=complex( slhaPtr->vmix(i,2), slhaPtr->imvmix(i,2) );
801 
802  // ~chi+ [i] ~chi- [j] Z : loop over [j]
803  for (int j=1; j<=nChar; j++) {
804 
805  // Chargino mixing
806  complex uj1, uj2, vj1, vj2;
807  uj1=complex( slhaPtr->umix(j,1), slhaPtr->imumix(j,1) );
808  uj2=complex( slhaPtr->umix(j,2), slhaPtr->imumix(j,2) );
809  vj1=complex( slhaPtr->vmix(j,1), slhaPtr->imvmix(j,1) );
810  vj2=complex( slhaPtr->vmix(j,2), slhaPtr->imvmix(j,2) );
811 
812  // ~chi+ [i] ~chi- [j] Z : couplings
813  OLp[i][j] = -vi1*conj(vj1) - 0.5*vi2*conj(vj2)
814  + ( (i == j) ? sin2W : 0.0);
815  ORp[i][j] = -conj(ui1)*uj1 - 0.5*conj(ui2)*uj2
816  + ( (i == j) ? sin2W : 0.0);
817 
818  if (DBSUSY) {
819  // tmp: verbose output
820  cout << " OL' [" << i << "][" << j << "] = "
821  << scientific << setw(10) << OLp[i][j]
822  << " OR' [" << i << "][" << j << "] = "
823  << scientific << setw(10) << ORp[i][j] << endl;
824  }
825  }
826 
827  // Loop over quark [l] flavour
828  for (int l=1;l<=3;l++) {
829 
830  // Set quark [l] masses
831  // Initial guess 0,0,0,mc,mb,mt with the latter from the PDT
832  double mul = particleDataPtr->m0(2*l);
833  double mdl = particleDataPtr->m0(2*l-1);
834  if (l == 1 || l == 2) { mul=0.0 ; mdl=0.0; }
835 
836  // Compute running mass from Yukawas and vevs if possible.
837  if (slhaPtr->yd.exists() && slhaPtr->hmix.exists(3)) {
838  double yll=slhaPtr->yd(l,l);
839  double v1=slhaPtr->hmix(3)/sqrt(1+pow(tanb,2));
840  if (yll > 0.0) mdl = yll * v1 / sqrt(2.0) ;
841  }
842  if (slhaPtr->yu.exists() && slhaPtr->hmix.exists(3)) {
843  double yll=slhaPtr->yu(l,l);
844  double v2=slhaPtr->hmix(3)/sqrt(1.0+1.0/pow(tanb,2));
845  if (yll > 0.0) mul = yll * v2 / sqrt(2.0) ;
846  }
847 
848  // Loop over squark [j] flavour
849  for (int j=1;j<=6;j++) {
850 
851  // Loop over off-diagonal quark [k] generation
852  for (int k=1;k<=3;k++) {
853 
854  // Set quark [k] masses
855  // Initial guess 0,0,0,0,mb,mt with the latter from the PDT
856  double muk = particleDataPtr->m0(2*k);
857  double mdk = particleDataPtr->m0(2*k-1);
858  if (k == 1) { muk=0.0 ; mdk=0.0; }
859  if (k == 2) { mdk=0.0 ; muk=0.0; }
860 
861  // Compute running mass from Yukawas and vevs if possible.
862  if (slhaPtr->yd.exists() && slhaPtr->hmix.exists(3)) {
863  double ykk=slhaPtr->yd(k,k);
864  double v1=slhaPtr->hmix(3)/sqrt(1+pow(tanb,2));
865  if (ykk > 0.0) mdk = ykk * v1 / sqrt(2.0) ;
866  }
867  if (slhaPtr->yu.exists() && slhaPtr->hmix.exists(3)) {
868  double ykk=slhaPtr->yu(k,k);
869  double v2=slhaPtr->hmix(3)/sqrt(1.0+1.0/pow(tanb,2));
870  if (ykk > 0.0) muk = ykk * v2 / sqrt(2.0) ;
871  }
872 
873  // CKM matrix (use Pythia one if no SLHA)
874  // (NB: could also try input one if no running one found, but
875  // would then need to compute from Wolfenstein)
876  complex Vlk=VCKMgen(l,k);
877  complex Vkl=VCKMgen(k,l);
878  if (slhaPtr->vckm.exists()) {
879  Vlk=complex(slhaPtr->vckm(l,k),slhaPtr->imvckm(l,k));
880  Vkl=complex(slhaPtr->vckm(k,l),slhaPtr->imvckm(k,l));
881  }
882 
883  // Squark mixing
884  complex Rdjk = complex(Rd(j,k), imRd(j,k) );
885  complex Rdjk3 = complex(Rd(j,k+3),imRd(j,k+3));
886  complex Rujk = complex(Ru(j,k), imRu(j,k) );
887  complex Rujk3 = complex(Ru(j,k+3),imRu(j,k+3));
888 
889 
890  // ~d[j] u[l] ~chi+[i]
891  LsduX[j][l][i] += (ui1*conj(Rdjk)
892  - mdk*ui2*conj(Rdjk3)/rt2/mW/cosb)*Vlk;
893  RsduX[j][l][i] -= mul*conj(vi2)*Vlk*conj(Rdjk)/rt2/mW/sinb;
894 
895  // ~u[j] d[l] ~chi+[i]
896  LsudX[j][l][i] += (conj(vi1)*Rujk
897  - muk*conj(vi2)*Rujk3/rt2/mW/sinb)*Vkl;
898  RsudX[j][l][i] -= mdl*ui2*Vkl*Rujk/rt2/mW/cosb;
899 
900  }
901 
902  if (DBSUSY) {
903  if (max(abs(LsduX[j][l][i]),abs(RsduX[j][l][i])) > 1e-6) {
904  // tmp: verbose output
905  cout << " LsduX[" << j << "][" << l << "][" << i << "] = "
906  << scientific << setw(10) << LsduX[j][l][i];
907  cout << " RsduX[" << j << "][" << l << "][" << i << "] = "
908  << scientific << setw(10) << RsduX[j][l][i] << endl;
909  }
910  if (max(abs(LsudX[j][l][i]),abs(RsudX[j][l][i])) > 1e-6) {
911  // tmp: verbose output
912  cout << " LsudX[" << j << "][" << l << "][" << i << "] = "
913  << scientific << setw(10) << LsudX[j][l][i];
914  cout << " RsudX[" << j << "][" << l << "][" << i << "] = "
915  << scientific << setw(10) << RsudX[j][l][i] << endl;
916  }
917  }
918  }
919  }
920 
921  // Loop over slepton [j] flavour
922  for (int j=1;j<=6;j++) {
923  for (int k=1;k<=3;k++) {
924 
925  LslvX[j][k][i] = 0.0;
926  RslvX[j][k][i] = 0.0;
927  LsvlX[j][k][i] = 0.0;
928  RsvlX[j][k][i] = 0.0;
929 
930  // Set lepton [k] masses
931  double ml = 0.0;
932  if (k == 3) ml = particleDataPtr->m0(15);
933 
934  if(j==k || j==k+3){ // No lepton mixing
935 
936  // ~l[j] v[l] ~chi+[i]
937  LslvX[j][k][i] += ui1- ml*ui2/rt2/mW/cosb;
938  RslvX[j][k][i] -= ml*conj(vi2)/rt2/mW/sinb;
939 
940  // ~v[j] l[l] ~chi+[i]
941  if(j<=3){ // No right handed sneutrinos
942  LsvlX[j][k][i] += conj(vi1) - ml*conj(vi2)/rt2/mW/sinb;
943  }
944  }
945 
946  if (DBSUSY) {
947  if (max(abs(LslvX[j][k][i]),abs(RslvX[j][k][i])) > 1e-6) {
948  // tmp: verbose output
949  cout << " LslvX[" << j << "][" << k << "][" << i << "] = "
950  << scientific << setw(10) << LslvX[j][k][i];
951  cout << " RslvX[" << j << "][" << k << "][" << i << "] = "
952  << scientific << setw(10) << RslvX[j][k][i] << endl;
953  }
954  if (max(abs(LsvlX[j][k][i]),abs(RsvlX[j][k][i])) > 1e-6) {
955  // tmp: verbose output
956  cout << " LsvlX[" << j << "][" << k << "][" << i << "] = "
957  << scientific << setw(10) << LsvlX[j][k][i];
958  cout << " RsvlX[" << j << "][" << k << "][" << i << "] = "
959  << scientific << setw(10) << RsvlX[j][k][i] << endl;
960  }
961  }
962  }
963  }
964  }
965 
966  // SLHA2 compatibility
967  // Check whether scalar particle masses are ordered
968  bool orderedQ = true;
969  bool orderedL = true;
970  for (int j=1;j<=5;j++) {
971  if (particleDataPtr->m0(idSlep(j+1)) < particleDataPtr->m0(idSlep(j)))
972  orderedL = false;
973  if (particleDataPtr->m0(idSup(j+1)) < particleDataPtr->m0(idSup(j)))
974  orderedQ = false;
975  if (particleDataPtr->m0(idSdown(j+1)) <particleDataPtr->m0(idSdown(j)))
976  orderedQ = false;
977  }
978 
979  // If ordered, change sparticle labels to mass-ordered enumeration
980  for (int i=1;i<=6;i++) {
981  ostringstream indx;
982  indx << i;
983  string uName = "~u_"+indx.str();
984  string dName = "~d_"+indx.str();
985  string lName = "~e_"+indx.str();
986  ParticleDataEntry* pdePtr;
987  if (orderedQ) {
988  pdePtr = particleDataPtr->particleDataEntryPtr(idSup(i));
989  pdePtr->setNames(uName,uName+"bar");
990  pdePtr = particleDataPtr->particleDataEntryPtr(idSdown(i));
991  pdePtr->setNames(dName,dName+"bar");
992  }
993  if (orderedL) {
994  pdePtr = particleDataPtr->particleDataEntryPtr(idSlep(i));
995  pdePtr->setNames(lName,lName+"bar");
996  }
997  }
998 
999  // Shorthand for RPV couplings
1000  // The input LNV lambda couplings
1001  LHtensor3Block<3> rvlle(slhaPtr->rvlamlle);
1002  // The input LNV lambda' couplings
1003  LHtensor3Block<3> rvlqd(slhaPtr->rvlamlqd);
1004  // The input BNV lambda'' couplings
1005  LHtensor3Block<3> rvudd(slhaPtr->rvlamudd);
1006 
1007  isLLE = false;
1008  isLQD = false;
1009  isUDD = false;
1010 
1011  // Symmetry properties
1012  if(rvlle.exists()){
1013  isLLE = true;
1014  for(int i=1;i<=3;i++){
1015  for(int j=i;j<=3;j++){
1016  //lambda(i,j,k)=-lambda(j,i,k)
1017  for(int k=1;k<=3;k++){
1018  if(i==j){
1019  rvLLE[i][j][k] = 0.0;
1020  }else{
1021  rvLLE[i][j][k] = rvlle(i,j,k);
1022  rvLLE[j][i][k] = -rvlle(i,j,k);
1023  }
1024  }
1025  }
1026  }
1027  }
1028 
1029  if(rvlqd.exists()){
1030  isLQD=true;
1031  for(int i=1;i<=3;i++){
1032  for(int j=1;j<=3;j++){
1033  for(int k=1;k<=3;k++){
1034  rvLQD[i][j][k] = rvlqd(i,j,k);
1035  }
1036  }
1037  }
1038  }
1039 
1040  //lambda''(k,j,i)=-lambda''(k,i,j)
1041 
1042  if(rvudd.exists()){
1043  isUDD = true;
1044  for(int i=1;i<=3;i++){
1045  for(int j=i;j<=3;j++){
1046  for(int k=1;k<=3;k++){
1047  if(i==j){
1048  rvUDD[k][i][j] = 0.0;
1049  }else{
1050  rvUDD[k][i][j] = rvudd(k,i,j);
1051  rvUDD[k][j][i] = -rvudd(k,i,j);
1052  }
1053  }
1054  }
1055  }
1056  }
1057 
1058  if(DBSUSY){
1059  for(int i=1;i<=3;i++){
1060  for(int j=1;j<=3;j++){
1061  for(int k=1;k<=3;k++){
1062  if(rvlle.exists())
1063  cout<<"LambdaLLE["<<i<<"]["<<j<<"]["<<k<<"]="<<rvLLE[i][j][k]<<" ";
1064  if(rvlqd.exists())
1065  cout<<"LambdaLQD["<<i<<"]["<<j<<"]["<<k<<"]="<<rvLQD[i][j][k]<<" ";
1066  if(rvudd.exists())
1067  cout<<"LambdaUDD["<<i<<"]["<<j<<"]["<<k<<"]="<<rvUDD[i][j][k]
1068  <<"\n";
1069  }
1070  }
1071  }
1072  }
1073 
1074  // Store the squark mixing matrix
1075  for(int i=1;i<=6;i++){
1076  for(int j=1;j<=3;j++){
1077  Rusq[i][j] = complex(Ru(i,j), imRu(i,j) );
1078  Rusq[i][j+3] = complex(Ru(i,j+3),imRu(i,j+3));
1079  Rdsq[i][j] = complex(Rd(i,j), imRd(i,j) );
1080  Rdsq[i][j+3] = complex(Rd(i,j+3),imRd(i,j+3));
1081  }
1082  }
1083 
1084  if(DBSUSY){
1085  cout<<"Ru"<<endl;
1086  for(int i=1;i<=6;i++){
1087  for(int j=1;j<=6;j++){
1088  cout << real(Rusq[i][j])<<" ";
1089  }
1090  cout<<endl;
1091  }
1092  cout<<"Rd"<<endl;
1093  for(int i=1;i<=6;i++){
1094  for(int j=1;j<=6;j++){
1095  cout << real(Rdsq[i][j])<<" ";
1096  }
1097  cout<<endl;
1098  }
1099  }
1100 
1101  // Let everyone know we are ready
1102  isInit = true;
1103 }
1104 
1105 //--------------------------------------------------------------------------
1106 
1107 // Return neutralino flavour codes.
1108 
1109 int CoupSUSY::idNeut(int idChi) {
1110 
1111  int id = 0;
1112  if (idChi == 1) id = 1000022;
1113  else if (idChi == 2) id = 1000023;
1114  else if (idChi == 3) id = 1000025;
1115  else if (idChi == 4) id = 1000035;
1116  else if (idChi == 5) id = 1000045;
1117  return id;
1118 }
1119 
1120 //--------------------------------------------------------------------------
1121 
1122 // Return chargino flavour codes.
1123 
1124 int CoupSUSY::idChar(int idChi) {
1125 
1126  int id = 0;
1127  if (idChi == 1) id = 1000024;
1128  else if (idChi == -1) id = -1000024;
1129  else if (idChi == 2) id = 1000037;
1130  else if (idChi == -2) id = -1000037;
1131  return id;
1132 }
1133 
1134 //--------------------------------------------------------------------------
1135 
1136 // Return sup flavour codes.
1137 
1138 int CoupSUSY::idSup(int iSup) {
1139 
1140  int id = 0;
1141  int sgn = ( iSup > 0 ) ? 1 : -1;
1142  iSup = abs(iSup);
1143  if (iSup == 1) id = 1000002;
1144  else if (iSup == 2) id = 1000004;
1145  else if (iSup == 3) id = 1000006;
1146  else if (iSup == 4) id = 2000002;
1147  else if (iSup == 5) id = 2000004;
1148  else if (iSup == 6) id = 2000006;
1149  return sgn*id;
1150 }
1151 
1152 //--------------------------------------------------------------------------
1153 
1154 // Return sdown flavour codes
1155 
1156 int CoupSUSY::idSdown(int iSdown) {
1157 
1158  int id = 0;
1159  int sgn = ( iSdown > 0 ) ? 1 : -1;
1160  iSdown = abs(iSdown);
1161  if (iSdown == 1) id = 1000001;
1162  else if (iSdown == 2) id = 1000003;
1163  else if (iSdown == 3) id = 1000005;
1164  else if (iSdown == 4) id = 2000001;
1165  else if (iSdown == 5) id = 2000003;
1166  else if (iSdown == 6) id = 2000005;
1167  return sgn*id;
1168 }
1169 
1170 //--------------------------------------------------------------------------
1171 
1172 // Function to return slepton flavour codes
1173 
1174 int CoupSUSY::idSlep(int iSlep) {
1175 
1176  int id = 0;
1177  int sgn = ( iSlep > 0 ) ? 1 : -1;
1178  iSlep = abs(iSlep);
1179  if (iSlep == 1) id = 1000011;
1180  else if (iSlep == 2) id = 1000013;
1181  else if (iSlep == 3) id = 1000015;
1182  else if (iSlep == 4) id = 2000011;
1183  else if (iSlep == 5) id = 2000013;
1184  else if (iSlep == 6) id = 2000015;
1185  return sgn*id;
1186 }
1187 
1188 //--------------------------------------------------------------------------
1189 
1190 // Return a particle name, given the PDG code.
1191 
1192 string CoupSUSY::getName(int pdgCode) {
1193 
1194  // Absolute value and corresponding SM code
1195  int codeA = abs(pdgCode);
1196  int idSM = codeA % 1000000;
1197 
1198  // Name
1199  string name;
1200 
1201  // Flag to indicate whether SLHA1 or SLHA2
1202  bool slha1 = false;
1203 
1204  // SM particles
1205  if (codeA == idSM) {
1206  // Neutrinos
1207  if (!slhaPtr->upmns.exists()) slha1=true;
1208  if (codeA == 12) name = (slha1) ? "nu_e" : "nu_1";
1209  if (codeA == 14) name = (slha1) ? "nu_mu" : "nu_2";
1210  if (codeA == 16) name = (slha1) ? "nu_tau" : "nu_3";
1211  }
1212 
1213  // Squarks
1214  else if ( idSM <= 6) {
1215  // up squarks
1216  if (idSM % 2 == 0) {
1217  // If SLHA1, return old PDG names
1218  if (slhaPtr->stopmix.exists()) slha1 = true;
1219  if (codeA == 1000002) name = (slha1) ? "~u_L" : "~u_1";
1220  if (codeA == 1000004) name = (slha1) ? "~c_L" : "~u_2";
1221  if (codeA == 1000006) name = (slha1) ? "~t_1" : "~u_3";
1222  if (codeA == 2000002) name = (slha1) ? "~u_R" : "~u_4";
1223  if (codeA == 2000004) name = (slha1) ? "~c_R" : "~u_5";
1224  if (codeA == 2000006) name = (slha1) ? "~t_2" : "~u_6";
1225  }
1226  // down squarks
1227  else {
1228  // If SLHA1, return old PDG names
1229  if (slhaPtr->sbotmix.exists()) slha1 = true;
1230  if (codeA == 1000001) name = (slha1) ? "~d_L" : "~d_1";
1231  if (codeA == 1000003) name = (slha1) ? "~s_L" : "~d_2";
1232  if (codeA == 1000005) name = (slha1) ? "~b_1" : "~d_3";
1233  if (codeA == 2000001) name = (slha1) ? "~d_R" : "~d_4";
1234  if (codeA == 2000003) name = (slha1) ? "~s_R" : "~d_5";
1235  if (codeA == 2000005) name = (slha1) ? "~b_2" : "~d_6";
1236  }
1237  if (pdgCode < 0) name += "bar";
1238  }
1239 
1240  // Sleptons
1241  else if ( idSM <= 19 ) {
1242 
1243  // Sneutrinos
1244  if (idSM % 2 == 0) {
1245  // If SLHA1, return old PDG names
1246  if (slhaPtr->staumix.exists()) slha1 = true;
1247  if (codeA == 1000012) name = (slha1) ? "~nu_eL" : "~nu_1";
1248  if (codeA == 1000014) name = (slha1) ? "~nu_muL" : "~nu_2";
1249  if (codeA == 1000016) name = (slha1) ? "~nu_tauL" : "~nu_3";
1250  if (codeA == 2000012) name = (slha1) ? "~nu_eR" : "~nu_4";
1251  if (codeA == 2000014) name = (slha1) ? "~nu_muR" : "~nu_5";
1252  if (codeA == 2000016) name = (slha1) ? "~nu_tauR" : "~nu_6";
1253  if (pdgCode < 0) name += "bar";
1254  }
1255  // charged sleptons
1256  else {
1257  // If SLHA1, return old PDG names
1258  if (slhaPtr->staumix.exists()) slha1 = true;
1259  if (codeA == 1000011) name = (slha1) ? "~e_L" : "~e_1";
1260  if (codeA == 1000013) name = (slha1) ? "~mu_L" : "~e_2";
1261  if (codeA == 1000015) name = (slha1) ? "~tau_1" : "~e_3";
1262  if (codeA == 2000011) name = (slha1) ? "~e_R" : "~e_4";
1263  if (codeA == 2000013) name = (slha1) ? "~mu_R" : "~e_5";
1264  if (codeA == 2000015) name = (slha1) ? "~tau_2" : "~e_6";
1265  if (pdgCode < 0) name += "-";
1266  else name += "+";
1267  }
1268  }
1269 
1270  else if (codeA == 1000021) name = "~g";
1271  else if (codeA == 1000022) name = "~chi_10";
1272  else if (codeA == 1000023) name = "~chi_20";
1273  else if (codeA == 1000024) name = (pdgCode > 0) ? "~chi_1+" : "~chi_1-";
1274  else if (codeA == 1000025) name = "~chi_30";
1275  else if (codeA == 1000035) name = "~chi_40";
1276  else if (codeA == 1000037) name = (pdgCode > 0) ? "~chi_2+" : "~chi_2-";
1277 
1278  return name;
1279 
1280 }
1281 
1282 //--------------------------------------------------------------------------
1283 
1284 // Return neutralino code; zero if not a neutralino
1285 
1286 int CoupSUSY::typeNeut(int idPDG) {
1287  int type = 0;
1288  int idAbs = abs(idPDG);
1289  if(idAbs == 1000022) type = 1;
1290  else if(idAbs == 1000023) type = 2;
1291  else if(idAbs == 1000025) type = 3;
1292  else if(idAbs == 1000035) type = 4;
1293  else if(isNMSSM && idAbs == 1000045) type = 5;
1294  return type;
1295 
1296 }
1297 
1298 //--------------------------------------------------------------------------
1299 
1300 // Check whether particle is a Chargino
1301 
1302 int CoupSUSY::typeChar(int idPDG) {
1303  int type = 0;
1304  if(abs(idPDG) == 1000024) type = 1;
1305  else if (abs(idPDG) == 1000037)type = 2;
1306  return type;
1307 }
1308 
1309 //==========================================================================
1310 
1311 } // end namespace Pythia8
1312 
1313 
1314 
Generated by   doxygen 1.8.5