PYTHIA8/pythia8175/src/SusyCouplings.cxx

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