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