PYTHIA8/pythia8145/src/SusyCouplings.cxx

   1 // SusyCouplings.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2010 Peter Skands, Torbjorn Sjostrand.
   3 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
   4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
   5
   6 // Function definitions (not found in the header) for the
   7 // supersymmetric couplings class.
   8
   9 #include "SusyCouplings.h"
  10
  11 namespace Pythia8 {
  12
  13 //==========================================================================
  14
  15 // The CoupSUSY class.
  16
  17 //--------------------------------------------------------------------------
  18
  19 // Constants: could be changed here if desired, but normally should not.
  20 // These are of technical nature, as described for each.
  21
  22 // Allow verbose printout for debug purposes.
  23   const bool CoupSUSY::DEBUG = false;
  24
  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   // Is NMSSM switched on?
  39   isNMSSM = (slhaPtr->modsel(3) != 1 ? false : true);
  40   settingsPtr->flag("SLHA:NMSSM",isNMSSM);
  41   int nNeut = (isNMSSM ? 5 : 4);
  42   int nChar = 2;
  43
  44   // Initialize pole masses
  45   mZpole    = particleDataPtr->m0(23);
  46   wZpole    = particleDataPtr->mWidth(23);
  47   mWpole    = particleDataPtr->m0(24);
  48   wWpole    = particleDataPtr->mWidth(24);
  49
  50   // Running masses and weak mixing angle
  51   // (default to pole values if no running available)
  52   mW        = mWpole;
  53   mZ        = mZpole;
  54   sin2W     = sin2thetaW();
  55
  56   if (settingsPtr->mode("SUSY:sin2thetaWMode") == 3) {
  57     // Possibility to force on-shell sin2W definition, mostly intended for
  58     // cross-checks
  59     sin2W     = 1.0 - pow(mW/mZ,2);
  60   }else if (settingsPtr->mode("SUSY:sin2thetaWMode") == 2){
  61     // Possibility to use running sin2W definition, derived from gauge
  62   // couplings in running SLHA blocks (normally defined at SUSY scale)
  63     if(slhaPtr->gauge.exists(1) && slhaPtr->gauge.exists(2)
  64        && slhaPtr->hmix.exists(3)) {
  65       double gp=slhaPtr->gauge(1);
  66       double g =slhaPtr->gauge(2);
  67       double v =slhaPtr->hmix(3);
  68       mW      = g * v / 2.0;
  69       mZ      = sqrt(pow(gp,2)+pow(g,2)) * v / 2.0;
  70       double tan2W   = pow2(gp)/pow2(g);
  71       if (DEBUG) cout << " tan2W = " << tan2W << endl;
  72       sin2W   = pow2(gp)/(pow2(g)+pow2(gp));
  73     } else {
  74       cout<<" GAUGE block not found in SLHA; using sin(thetaW) at mZ"<<endl;
  75     }
  76   }
  77
  78   //Overload the SM value of sin2thetaW
  79   s2tW = sin2W;
  80
  81   sinW = sqrt(sin2W);
  82   cosW = sqrt(1.0-sin2W);
  83
  84   // Tan(beta)
  85   // By default, use the running one in HMIX (if not found, use MINPAR)
  86   tanb = slhaPtr->hmix.exists(2) ? slhaPtr->hmix(2) : slhaPtr->minpar(3);
  87   cosb = sqrt( 1.0 / (1.0 + tanb*tanb) );
  88   sinb = sqrt(max(0.0,1.0-cosb*cosb));
  89
  90   // tmp : verbose output
  91   if (DEBUG) {
  92     cout << " sin2W(Q) = " << sin2W << "  mW(Q) = " << mW
  93          << "  mZ(Q) = " << mZ << endl;
  94     cout << " vev(Q) = " << slhaPtr->hmix(3) << " tanb(Q) = " << tanb
  95          << endl;
  96     for (int i=1;i<=3;i++) {
  97       for (int j=1;j<=3;j++) {
  98         cout << " VCKM  [" << i << "][" << j << "] = "
  99              << scientific << setw(10) << VCKMgen(i,j) << endl;
 100       }
 101     }
 102   }
 103
 104   // Shorthand for squark mixing matrices
 105   SusyLesHouches::MatrixBlock<6> Ru(slhaPtr->usqmix);
 106   SusyLesHouches::MatrixBlock<6> Rd(slhaPtr->dsqmix);
 107   SusyLesHouches::MatrixBlock<6> imRu(slhaPtr->imusqmix);
 108   SusyLesHouches::MatrixBlock<6> imRd(slhaPtr->imdsqmix);
 109
 110   // Construct ~g couplings
 111   for (int i=1 ; i<=6 ; i++) {
 112     for (int j=1 ; j<=3 ; j++) {
 113       LsddG[i][j] = complex( Rd(i,j)  ,  imRd(i,j));
 114       RsddG[i][j] = complex(-Rd(i,j+3), -imRd(i,j+3));
 115       LsuuG[i][j] = complex( Ru(i,j)  ,  imRu(i,j));
 116       RsuuG[i][j] = complex(-Ru(i,j+3), -imRu(i,j+3));
 117
 118       if (DEBUG) {
 119         cout << " Lsddg  [" << i << "][" << j << "] = "
 120              << scientific << setw(10) << LsddG[i][j]
 121              << " RsddG  [" << i << "][" << j  << "] = "
 122              << scientific << setw(10) << RsddG[i][j] << endl;
 123         cout << " Lsuug  [" << i << "][" << j << "] = "
 124              << scientific << setw(10) << LsuuG[i][j]
 125              << " RsuuG  [" << i << "][" << j  << "] = "
 126              << scientific << setw(10) << RsuuG[i][j] << endl;
 127       }
 128     }
 129   }
 130
 131   // Construct qqZ couplings
 132   for (int i=1 ; i<=6 ; i++) {
 133
 134     // q[i] q[i] Z (def with extra factor 2 compared to [Okun])
 135     LqqZ[i] = af(i) - 2.0*ef(i)*sin2W ;
 136     RqqZ[i] =                  - 2.0*ef(i)*sin2W ;
 137
 138     // tmp: verbose output
 139     if (DEBUG) {
 140       cout << " LqqZ  [" << i << "][" << i << "] = "
 141            << scientific << setw(10) << LqqZ[i]
 142            << " RqqZ  [" << i << "][" << i  << "] = "
 143            << scientific << setw(10) << RqqZ[i] << endl;
 144     }
 145   }
 146
 147   // Construct ~q~qZ couplings
 148   for (int i=1 ; i<=6 ; i++) {
 149
 150     // Squarks can have off-diagonal couplings as well
 151     for (int j=1 ; j<=6 ; j++) {
 152
 153       // ~d[i] ~d[j] Z
 154       LsdsdZ[i][j] = 0.0;
 155       RsdsdZ[i][j] = 0.0;
 156       for (int k=1;k<=3;k++) {
 157         complex Rdik  = complex(Rd(i,k),  imRd(i,k)  );
 158         complex Rdjk  = complex(Rd(j,k),  imRd(j,k)  );
 159         complex Rdik3 = complex(Rd(i,k+3),imRd(i,k+3));
 160         complex Rdjk3 = complex(Rd(j,k+3),imRd(j,k+3));
 161         LsdsdZ[i][j] += LqqZ[1] * (Rdik*conj(Rdjk));
 162         RsdsdZ[i][j] += RqqZ[1] * (Rdik3*conj(Rdjk3));
 163       }
 164
 165       // ~u[i] ~u[j] Z
 166       LsusuZ[i][j] = 0.0;
 167       RsusuZ[i][j] = 0.0;
 168       for (int k=1;k<=3;k++) {
 169         complex Ruik  = complex(Ru(i,k)  ,imRu(i,k)  );
 170         complex Rujk  = complex(Ru(j,k)  ,imRu(j,k)  );
 171         complex Ruik3 = complex(Ru(i,k+3),imRu(i,k+3));
 172         complex Rujk3 = complex(Ru(j,k+3),imRu(j,k+3));
 173         LsusuZ[i][j] += LqqZ[2] * (Ruik*conj(Rujk));
 174         RsusuZ[i][j] += RqqZ[2] * (Ruik3*conj(Rujk3));
 175       }
 176
 177     // tmp: verbose output
 178       if (DEBUG) {
 179         if (max(abs(LsdsdZ[i][j]),abs(RsdsdZ[i][j])) > 1e-6) {
 180           cout << " LsdsdZ[" << i << "][" << j << "] = "
 181                << scientific << setw(10) << LsdsdZ[i][j]
 182                << " RsdsdZ[" << i << "][" << j << "] = "
 183                << scientific << setw(10) << RsdsdZ[i][j] << endl;
 184         }
 185         if (max(abs(LsusuZ[i][j]),abs(RsusuZ[i][j]))> 1e-6) {
 186           cout << " LsusuZ[" << i << "][" << j << "] = "
 187                << scientific << setw(10) << LsusuZ[i][j]
 188                << " RsusuZ[" << i << "][" << j << "] = "
 189                << scientific << setw(10) << RsusuZ[i][j] << endl;
 190         }
 191       }
 192
 193     }
 194
 195   }
 196
 197   // Construct udW couplings
 198
 199   // Loop over up [i] and down [j] quark generation
 200   for (int i=1;i<=3;i++) {
 201     for (int j=1;j<=3;j++) {
 202
 203       // CKM matrix (use Pythia one if no SLHA)
 204       // (NB: could also try input one if no running one found, but
 205       // would then need to compute from Wolfenstein)
 206       complex Vij=VCKMgen(i,j);
 207       if (slhaPtr->vckm.exists()) {
 208         Vij=complex(slhaPtr->vckm(i,j),slhaPtr->imvckm(i,j));
 209       }
 210
 211       // u[i] d[j] W
 212       LudW[i][j] = sqrt(2.0) * cosW * Vij;
 213       RudW[i][j] = 0.0;
 214
 215       // tmp: verbose output
 216       if (DEBUG) {
 217         cout << " LudW  [" << i << "][" << j << "] = "
 218              << scientific << setw(10) << LudW[i][j]
 219              << " RudW  [" << i << "][" << j << "] = "
 220              << scientific << setw(10) << RudW[i][j] << endl;
 221       }
 222     }
 223   }
 224
 225   // Construct ~u~dW couplings
 226
 227   // Loop over ~u[k] and ~d[l] flavours
 228   for (int k=1;k<=6;k++) {
 229     for (int l=1;l<=6;l++) {
 230
 231       LsusdW[k][l]=0.0;
 232       RsusdW[k][l]=0.0;
 233
 234       // Loop over u[i] and d[j] flavours
 235       for (int i=1;i<=3;i++) {
 236         for (int j=1;j<=3;j++) {
 237
 238           // CKM matrix (use Pythia one if no SLHA)
 239           // (NB: could also try input one if no running one found, but
 240           // would then need to compute from Wolfenstein)
 241           complex Vij=VCKMgen(i,j);
 242           if (slhaPtr->vckm.exists()) {
 243             Vij=complex(slhaPtr->vckm(i,j),slhaPtr->imvckm(i,j));
 244           }
 245
 246           // ~u[k] ~d[l] W (add one term for each quark flavour i,j)
 247           complex Ruki = complex(Ru(k,i),imRu(k,i));
 248           complex Rdlj = complex(Rd(l,j),imRd(l,j));
 249           LsusdW[k][l] += sqrt(2.0) * cosW * Vij * Ruki * conj(Rdlj);
 250           RsusdW[k][l] += 0.0;
 251
 252         }
 253       }
 254
 255       // tmp: verbose output
 256       if (DEBUG) {
 257         if (max(abs(LsusdW[k][l]),abs(RsusdW[k][l]))> 1e-6) {
 258           cout << " LsusdW[" << k << "][" << l << "] = "
 259                << scientific << setw(10) << LsusdW[k][l]
 260                << " RsusdW[" << k << "][" << l << "] = "
 261                << scientific << setw(10) << RsusdW[k][l] << endl;
 262         }
 263       }
 264
 265     }
 266   }
 267
 268   // Now we come to the ones with really many indices
 269
 270   // Construct ~chi0 couplings (allow for 5 neutralinos in NMSSM)
 271   for (int i=1;i<=nNeut;i++) {
 272
 273     // Ni1, Ni2, Ni3, Ni4, Ni5
 274     complex ni1,ni2,ni3,ni4,ni5;
 275     if (not isNMSSM) {
 276       ni1=complex( slhaPtr->nmix(i,1), slhaPtr->imnmix(i,1) );
 277       ni2=complex( slhaPtr->nmix(i,2), slhaPtr->imnmix(i,2) );
 278       ni3=complex( slhaPtr->nmix(i,3), slhaPtr->imnmix(i,3) );
 279       ni4=complex( slhaPtr->nmix(i,4), slhaPtr->imnmix(i,4) );
 280       ni5=complex( 0.0, 0.0);
 281     } else {
 282       ni1=complex( slhaPtr->nmnmix(i,1), slhaPtr->imnmnmix(i,1) );
 283       ni2=complex( slhaPtr->nmnmix(i,2), slhaPtr->imnmnmix(i,2) );
 284       ni3=complex( slhaPtr->nmnmix(i,3), slhaPtr->imnmnmix(i,3) );
 285       ni4=complex( slhaPtr->nmnmix(i,4), slhaPtr->imnmnmix(i,4) );
 286       ni5=complex( slhaPtr->nmnmix(i,5), slhaPtr->imnmnmix(i,5) );
 287     }
 288
 289     // Change to positive mass convention
 290     complex iRot( 0., 1.);
 291     if (slhaPtr->mass(idNeut(i)) < 0.) {
 292       ni1 *= iRot;
 293       ni2 *= iRot;
 294       ni3 *= iRot;
 295       ni4 *= iRot;
 296       ni5 *= iRot;
 297     }
 298
 299     // ~chi0 [i] ~chi0 [j] Z : loop over [j]
 300     for (int j=1; j<=nNeut; j++) {
 301
 302       // neutralino [j] higgsino components
 303       complex nj3, nj4;
 304       if (not isNMSSM) {
 305         nj3=complex( slhaPtr->nmix(j,3), slhaPtr->imnmix(j,3) );
 306         nj4=complex( slhaPtr->nmix(j,4), slhaPtr->imnmix(j,4) );
 307       } else {
 308         nj3=complex( slhaPtr->nmnmix(j,3), slhaPtr->imnmnmix(j,3) );
 309         nj4=complex( slhaPtr->nmnmix(j,4), slhaPtr->imnmnmix(j,4) );
 310       }
 311       // Change to positive mass convention
 312       if (slhaPtr->mass(idNeut(j)) < 0.) {
 313         nj3 *= iRot;
 314         nj4 *= iRot;
 315       }
 316
 317       // ~chi0 [i] ~chi0 [j] Z : couplings
 318       OLpp[i][j] = -0.5 * ni3 * conj(nj3) + 0.5 * ni4 * conj(nj4);
 319       ORpp[i][j] = 0.5 * conj(ni3) * nj3 - 0.5 * conj(ni4) * nj4;
 320
 321     // tmp: verbose output
 322       if (DEBUG) {
 323         cout << " OL''  [" << i << "][" << j << "] = "
 324              << scientific << setw(10) << OLpp[i][j]
 325              << " OR''  [" << i << "][" << j << "] = "
 326              << scientific << setw(10) << ORpp[i][j] << endl;
 327       }
 328
 329     }
 330
 331     // ~chi0 [i] ~chi+ [j] W : loop over [j]
 332     for (int j=1; j<=nChar; j++) {
 333
 334       // Chargino mixing
 335       complex uj1, uj2, vj1, vj2;
 336       uj1=complex( slhaPtr->umix(j,1), slhaPtr->imumix(j,1) );
 337       uj2=complex( slhaPtr->umix(j,2), slhaPtr->imumix(j,2) );
 338       vj1=complex( slhaPtr->vmix(j,1), slhaPtr->imvmix(j,1) );
 339       vj2=complex( slhaPtr->vmix(j,2), slhaPtr->imvmix(j,2) );
 340
 341       // ~chi0 [i] ~chi+ [j] W : couplings
 342       OL[i][j] = -1.0/sqrt(2.0)*ni4*conj(vj2)+ni2*conj(vj1);
 343       OR[i][j] = 1.0/sqrt(2.0)*conj(ni3)*uj2+conj(ni2)*uj1;
 344
 345       // tmp: verbose output
 346       if (DEBUG) {
 347         cout << " OL    [" << i << "][" << j << "] = "
 348              << scientific << setw(10) << OL[i][j]
 349              << " OR    [" << i << "][" << j << "] = "
 350              << scientific << setw(10) << OR[i][j] << endl;
 351       }
 352     }
 353
 354     // Charges
 355     double ed  = -1.0/3.0;
 356     double T3d = -0.5;
 357     double eu  =  2.0/3.0;
 358     double T3u =  0.5;
 359
 360     // Loop over quark [k] generation
 361     for (int k=1;k<=3;k++) {
 362
 363       // Set quark masses
 364       // Initial guess 0,0,0,mc,mb,mt with the latter from the PDT
 365       double mu = particleDataPtr->m0(2*k);
 366       double md = particleDataPtr->m0(2*k-1);
 367       if (k == 1) { mu=0.0 ; md=0.0; }
 368       if (k == 2) { md=0.0 ; mu=0.0; }
 369
 370       // Compute running mass from Yukawas and vevs if possible.
 371       if (slhaPtr->yd.exists() && slhaPtr->hmix.exists(3)) {
 372         double ykk=slhaPtr->yd(k,k);
 373         double v1=slhaPtr->hmix(3)/sqrt(1+pow(tanb,2));
 374         if (ykk > 0.0) md = ykk * v1 / sqrt(2.0) ;
 375       }
 376       if (slhaPtr->yu.exists() && slhaPtr->hmix.exists(3)) {
 377         double ykk=slhaPtr->yu(k,k);
 378         double v2=slhaPtr->hmix(3)/sqrt(1.0+1.0/pow(tanb,2));
 379         if (ykk > 0.0) mu = ykk * v2 / sqrt(2.0) ;
 380       }
 381
 382       // tmp: verbose output
 383       if (DEBUG) {
 384         cout  <<  " Gen = " << k << " mu = " << mu << " md = " << md
 385               << " yUU,DD = " << slhaPtr->yu(k,k) << ","
 386               << slhaPtr->yd(k,k) << endl;
 387       }
 388
 389       // Loop over squark [j] flavour
 390       for (int j=1;j<=6;j++) {
 391
 392         // Squark mixing
 393         complex Rdjk  = complex(Rd(j,k),  imRd(j,k)  );
 394         complex Rdjk3 = complex(Rd(j,k+3),imRd(j,k+3));
 395         complex Rujk  = complex(Ru(j,k),  imRu(j,k)  );
 396         complex Rujk3 = complex(Ru(j,k+3),imRu(j,k+3));
 397
 398         // ~d[j] d[k] ~chi0[i]
 399         // Changed according to new notation
 400         LsddX[j][k][i] = ((ed-T3d)*sinW*ni1 + T3d*cosW*ni2)*conj(Rdjk)
 401           + md*cosW*ni3*conj(Rdjk3)/2.0/mW/cosb;
 402         RsddX[j][k][i] = -ed*sinW*conj(ni1)*conj(Rdjk3)
 403           + md*cosW*conj(ni3)*conj(Rdjk)/2.0/mW/cosb;
 404
 405         // ~u[j] u[k] ~chi0[i]
 406         LsuuX[j][k][i] = ((eu-T3u)*sinW*ni1 + T3u*cosW*ni2)*conj(Rujk)
 407           + mu*cosW*ni4*conj(Rujk3)/2.0/mW/sinb;
 408         RsuuX[j][k][i] = -eu*sinW*conj(ni1)*conj(Rujk3)
 409           + mu*cosW*conj(ni4)*conj(Rujk)/2.0/mW/sinb;
 410
 411         if (DEBUG) {
 412           if (abs(LsddX[j][k][i]) > 1e-6) {
 413             // tmp: verbose output
 414             cout << " LsddX[" << j << "][" << k << "][" << i << "] = "
 415                  << scientific << setw(10) << LsddX[j][k][i] << endl;
 416           }
 417           if (abs(RsddX[j][k][i]) > 1e-6) {
 418             // tmp: verbose output
 419             cout << " RsddX[" << j << "][" << k << "][" << i << "] = "
 420                  << scientific << setw(10) << RsddX[j][k][i] << endl;
 421           }
 422           if (abs(LsuuX[j][k][i]) > 1e-6) {
 423             // tmp: verbose output
 424             cout << " LsuuX[" << j << "][" << k << "][" << i << "] = "
 425                  << scientific << setw(10) << LsuuX[j][k][i] << endl;
 426           }
 427           if (abs(RsuuX[j][k][i]) > 1e-6) {
 428             // tmp: verbose output
 429             cout << " RsuuX[" << j << "][" << k << "][" << i << "] = "
 430                  << scientific << setw(10) << RsuuX[j][k][i] << endl;
 431           }
 432         }
 433       }
 434
 435     }
 436
 437   }
 438
 439   // Construct ~chi+ couplings
 440   for (int i=1;i<=nChar;i++) {
 441
 442     // Ui1, Ui2, Vi1, Vi2
 443     complex ui1,ui2,vi1,vi2;
 444     ui1=complex( slhaPtr->umix(i,1), slhaPtr->imumix(i,1) );
 445     ui2=complex( slhaPtr->umix(i,2), slhaPtr->imumix(i,2) );
 446     vi1=complex( slhaPtr->vmix(i,1), slhaPtr->imvmix(i,1) );
 447     vi2=complex( slhaPtr->vmix(i,2), slhaPtr->imvmix(i,2) );
 448
 449     // ~chi+ [i] ~chi- [j] Z : loop over [j]
 450     for (int j=1; j<=nChar; j++) {
 451
 452       // Chargino mixing
 453       complex uj1, uj2, vj1, vj2;
 454       uj1=complex( slhaPtr->umix(j,1), slhaPtr->imumix(j,1) );
 455       uj2=complex( slhaPtr->umix(j,2), slhaPtr->imumix(j,2) );
 456       vj1=complex( slhaPtr->vmix(j,1), slhaPtr->imvmix(j,1) );
 457       vj2=complex( slhaPtr->vmix(j,2), slhaPtr->imvmix(j,2) );
 458
 459       // ~chi+ [i] ~chi- [j] Z : couplings
 460       OLp[i][j] = -vi1*conj(vj1) - 0.5*vi2*conj(vj2)
 461         + ( (i == j) ? sin2W : 0.0);
 462       ORp[i][j] = -conj(ui1)*uj1 - 0.5*conj(ui2)*uj2
 463         + ( (i == j) ? sin2W : 0.0);
 464
 465       if (DEBUG) {
 466         // tmp: verbose output
 467         cout << " OL'   [" << i << "][" << j << "] = "
 468              << scientific << setw(10) << OLp[i][j]
 469              << " OR'   [" << i << "][" << j << "] = "
 470              << scientific << setw(10) << ORp[i][j] << endl;
 471       }
 472     }
 473
 474     // Loop over quark [l] flavour
 475     for (int l=1;l<=3;l++) {
 476
 477       // Set quark [l] masses
 478       // Initial guess 0,0,0,mc,mb,mt with the latter from the PDT
 479       double mul = particleDataPtr->m0(2*l);
 480       double mdl = particleDataPtr->m0(2*l-1);
 481       if (l == 1 || l == 2) { mul=0.0 ; mdl=0.0; }
 482
 483       // Compute running mass from Yukawas and vevs if possible.
 484       if (slhaPtr->yd.exists() && slhaPtr->hmix.exists(3)) {
 485         double yll=slhaPtr->yd(l,l);
 486         double v1=slhaPtr->hmix(3)/sqrt(1+pow(tanb,2));
 487         if (yll > 0.0) mdl = yll * v1 / sqrt(2.0) ;
 488       }
 489       if (slhaPtr->yu.exists() && slhaPtr->hmix.exists(3)) {
 490         double yll=slhaPtr->yu(l,l);
 491         double v2=slhaPtr->hmix(3)/sqrt(1.0+1.0/pow(tanb,2));
 492         if (yll > 0.0) mul = yll * v2 / sqrt(2.0) ;
 493       }
 494
 495       // Loop over squark [j] flavour
 496       for (int j=1;j<=6;j++) {
 497
 498         // Loop over off-diagonal quark [k] generation
 499         for (int k=1;k<=3;k++) {
 500
 501           // Set quark [k] masses
 502           // Initial guess 0,0,0,0,mb,mt with the latter from the PDT
 503           double muk = particleDataPtr->m0(2*k);
 504           double mdk = particleDataPtr->m0(2*k-1);
 505           if (k == 1) { muk=0.0 ; mdk=0.0; }
 506           if (k == 2) { mdk=0.0 ; muk=0.0; }
 507
 508           // Compute running mass from Yukawas and vevs if possible.
 509           if (slhaPtr->yd.exists() && slhaPtr->hmix.exists(3)) {
 510             double ykk=slhaPtr->yd(k,k);
 511             double v1=slhaPtr->hmix(3)/sqrt(1+pow(tanb,2));
 512             if (ykk > 0.0) mdk = ykk * v1 / sqrt(2.0) ;
 513           }
 514           if (slhaPtr->yu.exists() && slhaPtr->hmix.exists(3)) {
 515             double ykk=slhaPtr->yu(k,k);
 516             double v2=slhaPtr->hmix(3)/sqrt(1.0+1.0/pow(tanb,2));
 517             if (ykk > 0.0) muk = ykk * v2 / sqrt(2.0) ;
 518           }
 519
 520           // CKM matrix (use Pythia one if no SLHA)
 521           // (NB: could also try input one if no running one found, but
 522           // would then need to compute from Wolfenstein)
 523           complex Vlk=VCKMgen(l,k);
 524           complex Vkl=VCKMgen(k,l);
 525           if (slhaPtr->vckm.exists()) {
 526             Vlk=complex(slhaPtr->vckm(l,k),slhaPtr->imvckm(l,k));
 527             Vkl=complex(slhaPtr->vckm(k,l),slhaPtr->imvckm(k,l));
 528           }
 529
 530           // Squark mixing
 531           complex Rdjk  = complex(Rd(j,k),  imRd(j,k)  );
 532           complex Rdjk3 = complex(Rd(j,k+3),imRd(j,k+3));
 533           complex Rujk  = complex(Ru(j,k),  imRu(j,k)  );
 534           complex Rujk3 = complex(Ru(j,k+3),imRu(j,k+3));
 535
 536           // sqrt(2)
 537           double rt2 = sqrt(2.0);
 538
 539           // ~d[j] u[l] ~chi+[i]
 540           LsduX[j][l][i] += (ui1*conj(Rdjk)
 541                              - mdk*ui2*conj(Rdjk3)/rt2/mW/cosb)*Vlk;
 542           RsduX[j][l][i] -= mul*conj(vi2)*Vlk*conj(Rdjk)/rt2/mW/sinb;
 543
 544           // ~u[j] d[l] ~chi+[i]
 545           LsudX[j][l][i] += (conj(vi1)*Rujk
 546                              - muk*conj(vi2)*Rujk3/rt2/mW/sinb)*Vkl;
 547           RsudX[j][l][i] -= mdl*ui2*Vkl*Rujk/rt2/mW/cosb;
 548
 549         }
 550
 551         if (DEBUG) {
 552           if (max(abs(LsduX[j][l][i]),abs(RsduX[j][l][i])) > 1e-6) {
 553             // tmp: verbose output
 554             cout << " LsduX[" << j << "][" << l << "][" << i << "] = "
 555                  << scientific << setw(10) << LsduX[j][l][i];
 556             cout << " RsduX[" << j << "][" << l << "][" << i << "] = "
 557                  << scientific << setw(10) << RsduX[j][l][i] << endl;
 558           }
 559           if (max(abs(LsudX[j][l][i]),abs(RsudX[j][l][i])) > 1e-6) {
 560             // tmp: verbose output
 561             cout << " LsudX[" << j << "][" << l << "][" << i << "] = "
 562                  << scientific << setw(10) << LsudX[j][l][i];
 563             cout << " RsudX[" << j << "][" << l << "][" << i << "] = "
 564                  << scientific << setw(10) << RsudX[j][l][i] << endl;
 565           }
 566         }
 567
 568       }
 569
 570
 571     }
 572
 573   }
 574
 575   // SLHA2 compatibility
 576   // Check whether scalar particle masses are ordered
 577   bool orderedQ = true;
 578   bool orderedL = true;
 579   for (int j=1;j<=5;j++) {
 580     if (particleDataPtr->m0(idSlep(j+1)) < particleDataPtr->m0(idSlep(j)))
 581       orderedL  = false;
 582     if (particleDataPtr->m0(idSup(j+1)) < particleDataPtr->m0(idSup(j)))
 583       orderedQ  = false;
 584     if (particleDataPtr->m0(idSdown(j+1)) <particleDataPtr->m0(idSdown(j)))
 585       orderedQ  = false;
 586   }
 587
 588   // If ordered, change sparticle labels to mass-ordered enumeration
 589   for (int i=1;i<=6;i++) {
 590     ostringstream indx;
 591     indx << i;
 592     string uName = "~u_"+indx.str();
 593     string dName = "~d_"+indx.str();
 594     string lName = "~e_"+indx.str();
 595     ParticleDataEntry* pdePtr;
 596     if (orderedQ) {
 597       pdePtr = particleDataPtr->particleDataEntryPtr(idSup(i));
 598       pdePtr->setNames(uName,uName+"bar");
 599       pdePtr = particleDataPtr->particleDataEntryPtr(idSdown(i));
 600       pdePtr->setNames(dName,dName+"bar");
 601     }
 602     if (orderedL) {
 603       pdePtr = particleDataPtr->particleDataEntryPtr(idSlep(i));
 604       pdePtr->setNames(lName,lName+"bar");
 605     }
 606   }
 607
 608   // Shorthand for RPV couplings
 609   SusyLesHouches::Tensor3Block<3> rvlle(slhaPtr->rvlamlle); // The input LNV lambda couplings
 610   SusyLesHouches::Tensor3Block<3> rvlqd(slhaPtr->rvlamlqd); // The input LNV lambda' couplings
 611   SusyLesHouches::Tensor3Block<3> rvudd(slhaPtr->rvlamudd); // The input BNV lambda'' couplings
 612
 613   isLLE = false;
 614   isLQD = false;
 615   isUDD = false;
 616
 617   // Symmetry properties
 618   if(rvlle.exists()){
 619     isLLE = true;
 620     for(int i=1;i<=3;i++){
 621       for(int j=i;j<=3;j++){
 622         //lambda(i,j,k)=-lambda(j,i,k)
 623         for(int k=1;k<=3;k++){
 624           if(i==j){
 625             rvLLE[i][j][k] = 0.0;
 626           }else{
 627             rvLLE[i][j][k] = rvlle(i,j,k);
 628             rvLLE[j][i][k] = -rvlle(i,j,k);
 629           }
 630         }
 631       }
 632     }
 633   }
 634
 635   if(rvlqd.exists()){
 636     isLQD=true;
 637     for(int i=1;i<=3;i++){
 638       for(int j=1;j<=3;j++){
 639         for(int k=1;k<=3;k++){
 640             rvLQD[i][j][k] = rvlqd(i,j,k);
 641         }
 642       }
 643     }
 644   }
 645
 646   //lambda''(k,j,i)=-lambda''(k,i,j)
 647
 648   if(rvudd.exists()){
 649     isUDD = true;
 650     for(int i=1;i<=3;i++){
 651       for(int j=i;j<=3;j++){
 652         for(int k=1;k<=3;k++){
 653           if(i==j){
 654             rvUDD[k][i][j] = 0.0;
 655           }else{
 656             rvUDD[k][i][j] = rvudd(k,i,j);
 657             rvUDD[k][j][i] = -rvudd(k,i,j);
 658           }
 659         }
 660       }
 661     }
 662   }
 663
 664
 665   if(DEBUG){
 666     for(int i=1;i<=3;i++){
 667       for(int j=1;j<=3;j++){
 668         for(int k=1;k<=3;k++){
 669           if(rvlle.exists())
 670             cout<<"LambdaLLE["<<i<<"]["<<j<<"]["<<k<<"]="<<rvLLE[i][j][k]<<" ";
 671           if(rvlqd.exists())
 672             cout<<"LambdaLQD["<<i<<"]["<<j<<"]["<<k<<"]="<<rvLQD[i][j][k]<<" ";
 673           if(rvudd.exists())
 674             cout<<"LambdaUDD["<<i<<"]["<<j<<"]["<<k<<"]="<<rvUDD[i][j][k]<<"\n";
 675         }
 676       }
 677     }
 678   }
 679
 680   // Store the squark mixing matrix
 681   for(int i=1;i<=6;i++){
 682     for(int j=1;j<=3;j++){
 683       Rusq[i][j]   = complex(Ru(i,j),  imRu(i,j)  );
 684       Rusq[i][j+3] = complex(Ru(i,j+3),imRu(i,j+3));
 685       Rdsq[i][j]   = complex(Rd(i,j),  imRd(i,j)  );
 686       Rdsq[i][j+3] = complex(Rd(i,j+3),imRd(i,j+3));
 687     }
 688   }
 689
 690   if(DEBUG){
 691     cout<<"Ru"<<endl;
 692     for(int i=1;i<=6;i++){
 693       for(int j=1;j<=6;j++){
 694         cout << real(Rusq[i][j])<<"  ";
 695       }
 696       cout<<endl;
 697     }
 698     cout<<"Rd"<<endl;
 699     for(int i=1;i<=6;i++){
 700       for(int j=1;j<=6;j++){
 701         cout << real(Rdsq[i][j])<<"  ";
 702       }
 703       cout<<endl;
 704     }
 705   }
 706
 707   // Let everyone know we are ready
 708   isInit = true;
 709 }
 710
 711 //--------------------------------------------------------------------------
 712
 713 // Return neutralino flavour codes.
 714
 715 int CoupSUSY::idNeut(int idChi) {
 716   int id = 0;
 717   if      (idChi == 1) id = 1000022;
 718   else if (idChi == 2) id = 1000023;
 719   else if (idChi == 3) id = 1000025;
 720   else if (idChi == 4) id = 1000035;
 721   else if (idChi == 5) id = 1000045;
 722   return id;
 723 }
 724
 725 //--------------------------------------------------------------------------
 726
 727 // Return chargino flavour codes.
 728
 729 int CoupSUSY::idChar(int idChi) {
 730   int id = 0;
 731   if      (idChi ==  1) id =  1000024;
 732   else if (idChi == -1) id = -1000024;
 733   else if (idChi ==  2) id =  1000037;
 734   else if (idChi == -2) id = -1000037;
 735   return id;
 736 }
 737
 738 //--------------------------------------------------------------------------
 739
 740 // Return sup flavour codes.
 741
 742 int CoupSUSY::idSup(int iSup) {
 743   int id = 0;
 744   int sgn = ( iSup > 0 ) ? 1 : -1;
 745   iSup = abs(iSup);
 746   if      (iSup ==  1) id =  1000002;
 747   else if (iSup ==  2) id =  1000004;
 748   else if (iSup ==  3) id =  1000006;
 749   else if (iSup ==  4) id =  2000002;
 750   else if (iSup ==  5) id =  2000004;
 751   else if (iSup ==  6) id =  2000006;
 752   return sgn*id;
 753 }
 754
 755 //--------------------------------------------------------------------------
 756
 757 // Return sdown flavour codes
 758
 759 int CoupSUSY::idSdown(int iSdown) {
 760   int id = 0;
 761   int sgn = ( iSdown > 0 ) ? 1 : -1;
 762   iSdown = abs(iSdown);
 763   if      (iSdown ==  1) id =  1000001;
 764   else if (iSdown ==  2) id =  1000003;
 765   else if (iSdown ==  3) id =  1000005;
 766   else if (iSdown ==  4) id =  2000001;
 767   else if (iSdown ==  5) id =  2000003;
 768   else if (iSdown ==  6) id =  2000005;
 769   return sgn*id;
 770 }
 771
 772 //--------------------------------------------------------------------------
 773
 774 // Function to return slepton flavour codes
 775
 776 int CoupSUSY::idSlep(int iSlep) {
 777   int id = 0;
 778   int sgn = ( iSlep > 0 ) ? 1 : -1;
 779   iSlep = abs(iSlep);
 780   if      (iSlep ==  1) id =  1000011;
 781   else if (iSlep ==  2) id =  1000013;
 782   else if (iSlep ==  3) id =  1000015;
 783   else if (iSlep ==  4) id =  2000011;
 784   else if (iSlep ==  5) id =  2000013;
 785   else if (iSlep ==  6) id =  2000015;
 786   return sgn*id;
 787 }
 788
 789 //--------------------------------------------------------------------------
 790
 791 // Return a particle name, given the PDG code.
 792
 793 string CoupSUSY::getName(int pdgCode) {
 794
 795   // Absolute value and corresponding SM code
 796   int codeA = abs(pdgCode);
 797   int idSM  = codeA % 1000000;
 798
 799   // Name
 800   string name;
 801
 802   // Flag to indicate whether SLHA1 or SLHA2
 803   bool slha1 = false;
 804
 805   // SM particles
 806   if (codeA == idSM) {
 807     // Neutrinos
 808     if (!slhaPtr->upmns.exists()) slha1=true;
 809     if (codeA == 12) name = (slha1) ? "nu_e" : "nu_1";
 810     if (codeA == 14) name = (slha1) ? "nu_mu" : "nu_2";
 811     if (codeA == 16) name = (slha1) ? "nu_tau" : "nu_3";
 812   }
 813
 814   // Squarks
 815   else if ( idSM <= 6) {
 816     // up squarks
 817     if (idSM % 2 == 0) {
 818       // If SLHA1, return old PDG names
 819       if (slhaPtr->stopmix.exists()) slha1 = true;
 820       if (codeA == 1000002) name = (slha1) ? "~u_L" : "~u_1";
 821       if (codeA == 1000004) name = (slha1) ? "~c_L" : "~u_2";
 822       if (codeA == 1000006) name = (slha1) ? "~t_1" : "~u_3";
 823       if (codeA == 2000002) name = (slha1) ? "~u_R" : "~u_4";
 824       if (codeA == 2000004) name = (slha1) ? "~c_R" : "~u_5";
 825       if (codeA == 2000006) name = (slha1) ? "~t_2" : "~u_6";
 826     }
 827     // down squarks
 828     else {
 829       // If SLHA1, return old PDG names
 830       if (slhaPtr->sbotmix.exists()) slha1 = true;
 831       if (codeA == 1000001) name = (slha1) ? "~d_L" : "~d_1";
 832       if (codeA == 1000003) name = (slha1) ? "~s_L" : "~d_2";
 833       if (codeA == 1000005) name = (slha1) ? "~b_1" : "~d_3";
 834       if (codeA == 2000001) name = (slha1) ? "~d_R" : "~d_4";
 835       if (codeA == 2000003) name = (slha1) ? "~s_R" : "~d_5";
 836       if (codeA == 2000005) name = (slha1) ? "~b_2" : "~d_6";
 837     }
 838     if (pdgCode < 0) name += "bar";
 839   }
 840
 841   // Sleptons
 842   else if ( idSM <= 19 ) {
 843
 844     // Sneutrinos
 845    if (idSM % 2 == 0) {
 846       // If SLHA1, return old PDG names
 847       if (slhaPtr->staumix.exists()) slha1 = true;
 848       if (codeA == 1000012) name = (slha1) ? "~nu_eL" : "~nu_1";
 849       if (codeA == 1000014) name = (slha1) ? "~nu_muL" : "~nu_2";
 850       if (codeA == 1000016) name = (slha1) ? "~nu_tauL" : "~nu_3";
 851       if (codeA == 2000012) name = (slha1) ? "~nu_eR" : "~nu_4";
 852       if (codeA == 2000014) name = (slha1) ? "~nu_muR" : "~nu_5";
 853       if (codeA == 2000016) name = (slha1) ? "~nu_tauR" : "~nu_6";
 854       if (pdgCode < 0) name += "bar";
 855     }
 856     // charged sleptons
 857     else {
 858       // If SLHA1, return old PDG names
 859       if (slhaPtr->staumix.exists()) slha1 = true;
 860       if (codeA == 1000011) name = (slha1) ? "~e_L" : "~e_1";
 861       if (codeA == 1000013) name = (slha1) ? "~mu_L" : "~e_2";
 862       if (codeA == 1000015) name = (slha1) ? "~tau_1" : "~e_3";
 863       if (codeA == 2000011) name = (slha1) ? "~e_R" : "~e_4";
 864       if (codeA == 2000013) name = (slha1) ? "~mu_R" : "~e_5";
 865       if (codeA == 2000015) name = (slha1) ? "~tau_2" : "~e_6";
 866       if (pdgCode < 0) name += "-";
 867       else name += "+";
 868     }
 869   }
 870
 871   else if (codeA == 1000021) name = "~g";
 872   else if (codeA == 1000022) name = "~chi_10";
 873   else if (codeA == 1000023) name = "~chi_20";
 874   else if (codeA == 1000024) name = (pdgCode > 0) ? "~chi_1+" : "~chi_1-";
 875   else if (codeA == 1000025) name = "~chi_30";
 876   else if (codeA == 1000035) name = "~chi_40";
 877   else if (codeA == 1000037) name = (pdgCode > 0) ? "~chi_2+" : "~chi_2-";
 878
 879   return name;
 880
 881 }
 882
 883
 884 //==========================================================================
 885
 886 } // end namespace Pythia8
 887
 888