PYTHIA8/pythia8145/src/SusyResonanceWidths.cxx

   1 // ResonanceWidths.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2010
   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
   7 // the ResonanceWidths class and classes derived from it.
   8
   9 #include "SusyResonanceWidths.h"
  10 #include "SusyCouplings.h"
  11 #include "PythiaComplex.h"
  12
  13 namespace Pythia8 {
  14
  15 //==========================================================================
  16
  17 // The SUSYResonanceWidths Class
  18 // Derived class for SUSY resonances
  19
  20 const bool SUSYResonanceWidths::DEBUG = false;
  21
  22 //--------------------------------------------------------------------------
  23
  24 bool SUSYResonanceWidths::init(Info* infoPtrIn, Settings* settingsPtrIn,
  25    ParticleData* particleDataPtrIn, Couplings* couplingsPtrIn) {
  26
  27   // Save pointers.
  28   infoPtr         = infoPtrIn;
  29   settingsPtr     = settingsPtrIn;
  30   particleDataPtr = particleDataPtrIn;
  31   coupSUSYPtr     = (couplingsPtrIn->isSUSY ? (CoupSUSY *) couplingsPtrIn: 0 );
  32
  33   // No initialization necessary for SM-only
  34   if(!couplingsPtrIn->isSUSY) return true;
  35
  36   // Minimal decaying-resonance width. Minimal phase space for meMode = 103.
  37   minWidth     = settingsPtr->parm("ResonanceWidths:minWidth");
  38   minThreshold = settingsPtr->parm("ResonanceWidths:minThreshold");
  39
  40   // Pointer to particle species.
  41   particlePtr  = particleDataPtr->particleDataEntryPtr(idRes);
  42   if (particlePtr == 0) {
  43     infoPtr->errorMsg("Error in ResonanceWidths::init:"
  44       " unknown resonance identity code");
  45     return false;
  46   }
  47
  48   // Generic particles should not have meMode < 100, but allow
  49   // some exceptions: not used Higgses and not used Technicolor.
  50   if (idRes == 35 || idRes == 36 || idRes == 37
  51     || idRes/1000000 == 3) isGeneric = false;
  52
  53   // Resonance properties: antiparticle, mass, width
  54   hasAntiRes   = particlePtr->hasAnti();
  55   mRes         = particlePtr->m0();
  56   GammaRes     = particlePtr->mWidth();
  57   m2Res        = mRes*mRes;
  58
  59   // For very narrow resonances assign fictitious small width.
  60   if (GammaRes < minWidth) GammaRes = 0.1 * minWidth;
  61   GamMRat      = GammaRes / mRes;
  62
  63   // Secondary decay chains by default all on.
  64   openPos      = 1.;
  65   openNeg      = 1.;
  66
  67   // Allow option where on-shell width is forced to current value.
  68   doForceWidth = particlePtr->doForceWidth();
  69   forceFactor  = 1.;
  70
  71   // Check that resonance OK.
  72   if (particlePtr == 0) infoPtr->errorMsg("Error in ResonanceWidths::init:"
  73       " unknown resonance identity code");
  74
  75   // Calculate various common prefactors for the current mass.
  76   mHat          = mRes;
  77
  78   // Initialize constants used for a resonance.
  79
  80   initConstants();
  81   calcPreFac(true);
  82
  83   // Reset quantities to sum. Declare variables inside loop.
  84   double widTot = 0.;
  85   double widPos = 0.;
  86   double widNeg = 0.;
  87   int    idNow, idAnti;
  88   double openSecPos, openSecNeg;
  89
  90   // Loop over all decay channels. Basic properties of channel.
  91   for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
  92     iChannel    = i;
  93     onMode      = particlePtr->channel(i).onMode();
  94     meMode      = particlePtr->channel(i).meMode();
  95     mult        = particlePtr->channel(i).multiplicity();
  96     widNow      = 0.;
  97
  98     // Warn if not relevant meMode.
  99     if ( meMode < 0 || meMode > 103 || (isGeneric && meMode < 100) ) {
 100       infoPtr->errorMsg("Error in ResonanceWidths::init:"
 101         " resonance meMode not acceptable");
 102     }
 103
 104     // Check if decay table was read in via SLHA
 105     bool hasDecayTable = false;
 106     for(unsigned int iDec = 1; iDec < (coupSUSYPtr->slhaPtr)->decays.size(); iDec++)
 107       hasDecayTable = ((coupSUSYPtr->slhaPtr)->decays[iDec].getId() == abs(idRes));
 108
 109     // Calculation of SUSY particle widths
 110     if (meMode == 103 && GammaRes > 0. &&
 111         (!settingsPtr->flag("SLHA:useDecayTable") || !hasDecayTable)) {
 112       // Read out information on channel: primarily use first two.
 113       id1       = particlePtr->channel(i).product(0);
 114       id2       = particlePtr->channel(i).product(1);
 115       id1Abs    = abs(id1);
 116       id2Abs    = abs(id2);
 117
 118       // Order first two in descending order of absolute values.
 119       if (id2Abs > id1Abs) {swap( id1, id2); swap( id1Abs, id2Abs);}
 120
 121       // Allow for third product to be treated in derived classes.
 122       if (mult > 2) {
 123         id3     = particlePtr->channel(i).product(2);
 124         id3Abs  = abs(id3);
 125
 126         // Also order third into descending order of absolute values.
 127         if (id3Abs > id2Abs) {swap( id2, id3); swap( id2Abs, id3Abs);}
 128         if (id2Abs > id1Abs) {swap( id1, id2); swap( id1Abs, id2Abs);}
 129       }
 130
 131       // Read out masses. Calculate two-body phase space.
 132       mf1       = particleDataPtr->m0(id1Abs);
 133       mf2       = particleDataPtr->m0(id2Abs);
 134       mr1       = pow2(mf1 / mHat);
 135       mr2       = pow2(mf2 / mHat);
 136       ps        = (mHat < mf1 + mf2 + MASSMARGIN) ? 0.
 137                 : sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 );
 138       if (mult > 2) {
 139         mf3     = particleDataPtr->m0(id3Abs);
 140         mr3     = pow2(mf3 / mHat);
 141       }
 142
 143       // Let derived class calculate width for channel provided.
 144       calcWidth(true);
 145     }
 146
 147     // Channels with meMode >= 100 are calculated based on stored values.
 148     else widNow = GammaRes * particlePtr->channel(i).bRatio();
 149
 150     // Find secondary open fractions of partial width.
 151     openSecPos  = 1.;
 152     openSecNeg  = 1.;
 153     if (widNow > 0.) for (int j = 0; j < mult; ++j) {
 154       idNow     = particlePtr->channel(i).product(j);
 155       idAnti    = (particleDataPtr->hasAnti(idNow)) ? -idNow : idNow;
 156       // Secondary widths not yet initialized for heavier states,
 157       // so have to assume unit open fraction there.
 158       if (idNow == 23 || abs(idNow) == 24
 159         || particleDataPtr->m0(abs(idNow)) < mRes) {
 160         openSecPos *= particleDataPtr->resOpenFrac(idNow);
 161         openSecNeg *= particleDataPtr->resOpenFrac(idAnti);
 162       }
 163     }
 164
 165     // Store partial widths and secondary open fractions.
 166     particlePtr->channel(i).onShellWidth(widNow);
 167     particlePtr->channel(i).openSec( idRes, openSecPos);
 168     particlePtr->channel(i).openSec(-idRes, openSecNeg);
 169
 170     // Update sum over all channnels and over open channels only.
 171     widTot     += widNow;
 172     if (onMode == 1 || onMode == 2) widPos += widNow * openSecPos;
 173     if (onMode == 1 || onMode == 3) widNeg += widNow * openSecNeg;
 174   }
 175
 176   // If no decay channels are open then set particle stable and done.
 177   if (widTot < minWidth) {
 178     particlePtr->setMayDecay(false, false);
 179     particlePtr->setMWidth(0., false);
 180     for (int i = 0; i < particlePtr->sizeChannels(); ++i)
 181       particlePtr->channel(i).bRatio( 0., false);
 182     return true;
 183   }
 184
 185   // Normalize branching ratios to unity.
 186   double bRatio;
 187   for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
 188     bRatio      = particlePtr->channel(i).onShellWidth() / widTot;
 189     particlePtr->channel(i).bRatio( bRatio, false);
 190   }
 191
 192   // Optionally force total width by rescaling of all partial ones.
 193   if (doForceWidth) {
 194     forceFactor = GammaRes / widTot;
 195     for (int i = 0; i < particlePtr->sizeChannels(); ++i)
 196       particlePtr->channel(i).onShellWidthFactor( forceFactor);
 197   }
 198
 199   // Else update newly calculated partial width.
 200   else {
 201     particlePtr->setMWidth(widTot, false);
 202     GammaRes    = widTot;
 203   }
 204
 205   // Updated width-to-mass ratio. Secondary widths for open.
 206   GamMRat       = GammaRes / mRes;
 207   openPos       = widPos / widTot;
 208   openNeg       = widNeg / widTot;
 209
 210   // Done.
 211   return true;
 212
 213 }
 214
 215
 216
 217 //==========================================================================
 218
 219 // The ResonanceSquark class
 220 // Derived class for Squark resonances
 221
 222 //--------------------------------------------------------------------------
 223
 224 // Initialize constants.
 225
 226 void ResonanceSquark::initConstants() {
 227
 228   // Locally stored properties and couplings.
 229   alpS  = coupSUSYPtr->alphaS(mHat * mHat );
 230   alpEM = coupSUSYPtr->alphaEM(mHat * mHat);
 231   s2W   = coupSUSYPtr->sin2W;
 232 }
 233
 234 //--------------------------------------------------------------------------
 235
 236 // Calculate various common prefactors for the current mass.
 237
 238 void ResonanceSquark::calcPreFac(bool) {
 239
 240   // Common coupling factors.
 241   preFac = 1.0 / (s2W * pow(mHat,3));
 242
 243 }
 244
 245 //--------------------------------------------------------------------------
 246
 247 // Calculate width for currently considered channel.
 248
 249 void ResonanceSquark::calcWidth(bool) {
 250
 251   // Squark type -- in u_i/d_i and generation
 252   int ksusy = 1000000;
 253   bool idown = (abs(idRes)%2 == 0 ? false : true);
 254   int isq = (abs(idRes)/ksusy == 2) ?
 255     (abs(idRes)%10+1)/2 + 3: (abs(idRes)%10+1)/2;
 256   int isqgen = (abs(idRes)%10 + 1)/2;
 257
 258   // Check that above threshold.
 259   if (ps == 0.) return;
 260   kinFac = (mHat * mHat - mf1 * mf1 -mf2 * mf2);
 261   lambda = lam(mHat*mHat, mf1*mf1, mf2*mf2);
 262
 263   double fac = 0.0 , wid = 0.0;
 264
 265   //Case 1: RPV decay
 266   if(id1Abs < 7 && id2Abs < 7){
 267
 268     //Temporary till colour assignments sorted out
 269     //widNow = 0.0;
 270     //return;
 271
 272     // Quark generations
 273     int iq1 = (id1Abs + 1)/2;
 274     int iq2 = (id2Abs + 1)/2;
 275
 276     // Check for RPV UDD couplings
 277     if(!coupSUSYPtr->isUDD) {widNow = 0; return;}
 278
 279     // ~q -> q_i + q_j
 280
 281     fac = mHat;
 282     if(idown) {
 283       if ((id1Abs+id2Abs)%2 == 1){
 284         if(id1Abs%2==1)
 285           for(int isq2 = 4; isq2 < 7; isq2++)
 286             wid = coupSUSYPtr->rvUDD[iq2][iq1][isqgen] *
 287               norm(coupSUSYPtr->Rdsq[isq][isq2]);
 288         else
 289           for(int isq2 = 4; isq2 < 7; isq2++)
 290             wid = coupSUSYPtr->rvUDD[iq1][iq2][isqgen] *
 291               norm(coupSUSYPtr->Rdsq[isq][isq2]);
 292       }
 293     }
 294     else {
 295       if ((id1Abs+id2Abs)%2 != 0) widNow = 0.0;
 296       else
 297         for(int isq2 = 4; isq2 < 7; isq2++)
 298           wid = coupSUSYPtr->rvUDD[isq][iq1][iq2] *
 299             norm(coupSUSYPtr->Rusq[isq][isq2]);
 300     }
 301   }
 302
 303   //Case 2: quark + gaugino (higgsino)
 304   else if (id1Abs > ksusy && id2Abs < 7) {
 305
 306     int iq = (id2Abs + 1)/2;
 307
 308     // ~q -> ~g + q
 309     if(id1Abs == 1000021 && idRes%10 == id2Abs) {
 310       fac = 2.0 / 3.0 * alpS *  preFac * sqrt(lambda);
 311       if(idown)
 312         wid = kinFac * (norm(coupSUSYPtr->LsddG[isq][iq])+norm(coupSUSYPtr->RsddG[isq][iq]))
 313             - 4.0 * mf1 * mf2 * real(coupSUSYPtr->LsddG[isq][iq]*coupSUSYPtr->RsddG[isq][iq]);
 314       else
 315         wid = kinFac * (norm(coupSUSYPtr->LsuuG[isq][iq])+norm(coupSUSYPtr->RsuuG[isq][iq]))
 316             - 4.0 * mf1 * mf2 * real(coupSUSYPtr->LsuuG[isq][iq]*coupSUSYPtr->RsuuG[isq][iq]);
 317
 318     }
 319     else
 320       for(int i=1; i<6 ; i++){
 321         // ~q -> ~chi0 + q
 322         if(coupSUSYPtr->idNeut(i)==id1Abs && idRes%2 == id2Abs%2){
 323           fac = alpEM *  preFac * sqrt(lambda)/ (4.0 * (1 - s2W));
 324           if(idown)
 325             wid = kinFac * (norm(coupSUSYPtr->LsddX[isq][iq][i]) + norm(coupSUSYPtr->RsddX[isq][iq][i]))
 326                 - 4.0 * mf1 * mf2 * real(coupSUSYPtr->LsddX[isq][iq][i]*coupSUSYPtr->RsddX[isq][iq][i]);
 327           else
 328             wid = kinFac * (norm(coupSUSYPtr->LsuuX[isq][iq][i]) + norm(coupSUSYPtr->RsuuX[isq][iq][i]))
 329                 - 4.0 * mf1 * mf2 * real(coupSUSYPtr->LsuuX[isq][iq][i]*coupSUSYPtr->RsuuX[isq][iq][i]);
 330         }
 331
 332         // ~q -> chi- + q
 333         else if (i < 3 && coupSUSYPtr->idChar(i)==id1Abs && idRes%2 != id2Abs%2){
 334
 335           fac = alpEM *  preFac * sqrt(lambda)/ (4.0 * (1 - s2W));
 336           if(idown)
 337             wid = kinFac * (norm(coupSUSYPtr->LsduX[isq][iq][i]) + norm(coupSUSYPtr->RsduX[isq][iq][i]))
 338                 - 4.0 * mf1 * mf2 * real(coupSUSYPtr->LsduX[isq][iq][i]*coupSUSYPtr->RsduX[isq][iq][i]);
 339           else
 340             wid = kinFac * (norm(coupSUSYPtr->LsudX[isq][iq][i]) + norm(coupSUSYPtr->RsudX[isq][iq][i]))
 341                 - 4.0 * mf1 * mf2 * real(coupSUSYPtr->LsudX[isq][iq][i]*coupSUSYPtr->RsudX[isq][iq][i]);
 342         }
 343       }
 344   }
 345
 346   //Case 3: ~q_i -> ~q_j + Z/W
 347   else if (id1Abs > ksusy && id1Abs%100 < 7 && (id2Abs == 23 || id2Abs == 24)){
 348
 349     fac = alpEM * preFac/(16.0 * pow2(particleDataPtr->m0(id2Abs)) * (1.0 - s2W))
 350       * pow(lambda, 1.5);
 351
 352     int isq2 = (id1Abs/ksusy == 2) ? (id1Abs%10+1)/2 + 3: (id1Abs%10+1)/2;
 353
 354     if(id2Abs == 23 && id1Abs%2 == idRes%2){
 355       if(idown)
 356         wid = norm(coupSUSYPtr->LsdsdZ[isq][isq2] + coupSUSYPtr->RsdsdZ[isq][isq2]);
 357       else
 358         wid = norm(coupSUSYPtr->LsusuZ[isq][isq2] + coupSUSYPtr->RsusuZ[isq][isq2]);
 359     }
 360     else if (id2Abs == 24 && id1Abs%2 != idRes%2){
 361       if(idown)
 362         wid = norm(coupSUSYPtr->LsusdW[isq2][isq]);
 363       else
 364         wid = norm(coupSUSYPtr->LsusdW[isq][isq2]);
 365     }
 366   }
 367
 368   widNow = fac * wid * ps;
 369   if(DEBUG) cout<<idRes<<":: id1:"<<id1Abs<<" id2:"<<id2Abs<<" Width: "<<widNow<<endl;
 370   return;
 371
 372 }
 373
 374 //==========================================================================
 375
 376 // The ResonanceGluino class
 377 // Derived class for Gluino resonances
 378
 379 //--------------------------------------------------------------------------
 380
 381 // Initialize constants.
 382
 383 void ResonanceGluino::initConstants() {
 384
 385   // Locally stored properties and couplings.
 386   alpS  = coupSUSYPtr->alphaS(mHat * mHat );
 387   return;
 388 }
 389
 390 //--------------------------------------------------------------------------
 391
 392 // Calculate various common prefactors for the current mass.
 393
 394 void ResonanceGluino::calcPreFac(bool) {
 395   // Common coupling factors.
 396   preFac = alpS /( 8.0 * pow(mHat,3));
 397   return;
 398 }
 399
 400 //--------------------------------------------------------------------------
 401
 402 // Calculate width for currently considered channel.
 403
 404 void ResonanceGluino::calcWidth(bool) {
 405
 406
 407   widNow = 0.0;
 408   if(ps == 0.) return;
 409   kinFac = (mHat * mHat - mf1 * mf1 + mf2 * mf2);
 410   lambda = lam(mHat*mHat, mf1*mf1 , mf2*mf2);
 411
 412   if(id1Abs > 1000000 && (id1Abs % 100) < 7 && id2Abs < 7) {
 413
 414     int isq = (abs(id1Abs)/1000000 == 2) ?
 415       (abs(id1Abs)%10+1)/2 + 3: (abs(id1Abs)%10+1)/2;
 416     bool idown = (id2Abs%2 == 1);
 417     int iq = (id2Abs + 1)/2;
 418
 419     // ~g -> ~q + q
 420     if(idown){
 421       widNow = kinFac * (norm(coupSUSYPtr->LsddG[isq][iq])+norm(coupSUSYPtr->RsddG[isq][iq]))
 422         + 4.0 * mHat * mf2 * real( coupSUSYPtr->LsddG[isq][iq]
 423                                                   * conj(coupSUSYPtr->RsddG[isq][iq]));
 424     }
 425     else{
 426       widNow = kinFac * (norm(coupSUSYPtr->LsuuG[isq][iq])+norm(coupSUSYPtr->RsuuG[isq][iq]))
 427         + 4.0 * mHat * mf2 * real( coupSUSYPtr->LsuuG[isq][iq]
 428                                                   * conj(coupSUSYPtr->RsuuG[isq][iq]));
 429
 430     }
 431     widNow = widNow * preFac * ps * sqrt(lambda);
 432     if(DEBUG) {
 433       cout<<"Gluino:: id1:"<<id1Abs<<" id2:"<<id2Abs<<" Width: ";
 434       cout<<scientific<<widNow<<endl;
 435     }
 436     return;
 437   }
 438 }
 439
 440 //==========================================================================
 441
 442 //  Class ResonanceNeut
 443 //  Derived class for Neutralino Resonances
 444
 445 //--------------------------------------------------------------------------
 446
 447
 448 void ResonanceNeut::initConstants(){
 449
 450   alpEM = coupSUSYPtr->alphaEM(mHat * mHat);
 451   s2W   = coupSUSYPtr->sin2W;
 452 }
 453
 454 //--------------------------------------------------------------------------
 455
 456 // Calculate various common prefactors for the current mass.
 457 void  ResonanceNeut::calcPreFac(bool){
 458
 459   // Common coupling factors.
 460   preFac = alpEM / (8.0 * s2W * pow(mHat,3));
 461   return;
 462 }
 463
 464 //--------------------------------------------------------------------------
 465
 466 // Calculate width for currently considered channel.
 467 void  ResonanceNeut::calcWidth(bool){
 468
 469   widNow = 0.0;
 470
 471   if(ps ==0.) return;
 472   kinFac = mHat * mHat - mf1 * mf1 + mf2 * mf2;
 473   kinFac2 = pow(mHat,4) + pow(mf1,4) - 2.0 * pow(mf2,4)
 474     + pow2(mHat) * pow2(mf2) + pow2(mf1) * pow2(mf2) - pow2(mHat) * pow2(mf1);
 475   lambda = lam(mHat*mHat, mf1*mf1, mf2*mf2);
 476
 477   // Stable lightest neutralino
 478   if(idRes == 1000022) return;
 479
 480   double fac = 0.0;
 481   int iNeut1 = typeNeut(idRes);
 482   int iNeut2 = typeNeut(id1Abs);
 483   int iChar1 = typeChar(id1Abs);
 484
 485   if(iNeut2>0 && id2Abs == 23){
 486     // ~chi0_i -> chi0_j + Z
 487     fac = kinFac2 * (norm(coupSUSYPtr->OLpp[iNeut1][iNeut2]) + norm(coupSUSYPtr->ORpp[iNeut1][iNeut2]));
 488     fac -= 12.0 * mHat * mf1 * pow2(mf2) *
 489       real(coupSUSYPtr->OLpp[iNeut1][iNeut2] * conj(coupSUSYPtr->ORpp[iNeut1][iNeut2]));
 490     fac /= pow2(mf2) * (1.0 - s2W);
 491   }
 492   else if(iChar1>0 && id2Abs==24){
 493     // ~chi0_i -> chi+_j + W- (or c.c.)
 494
 495     fac = kinFac2 * (norm(coupSUSYPtr->OL[iNeut1][iChar1]) + norm(coupSUSYPtr->OR[iNeut1][iChar1]));
 496     fac -= 12.0 * mHat * mf1 * pow2(mf2) *
 497       real(coupSUSYPtr->OL[iNeut1][iChar1] * conj(coupSUSYPtr->OR[iNeut1][iChar1]));
 498     fac /= pow2(mf2);
 499   }
 500   else if(id1Abs > 1000000 && id1Abs%100 < 7 && id2Abs < 7){
 501     // ~chi0_k -> ~q + q
 502     bool idown = (id1Abs%2 == 1);
 503     int iq = (id2Abs + 1 )/ 2;
 504     int isq = (abs(idRes)/1000000 == 2) ?
 505       (abs(idRes)%10+1)/2 + 3: (abs(idRes)%10+1)/2;
 506
 507     if(idown){
 508       fac = kinFac * (norm(coupSUSYPtr->LsddX[isq][iq][iNeut1]) + norm(coupSUSYPtr->RsddX[isq][iq][iNeut1]));
 509       fac += 4.0 * mHat * mf2 *
 510         real(coupSUSYPtr->LsddX[isq][iq][iNeut1] * conj(coupSUSYPtr->RsddX[isq][iq][iNeut1]));
 511     }
 512     else{
 513       fac = kinFac * (norm(coupSUSYPtr->LsuuX[isq][iq][iNeut1]) + norm(coupSUSYPtr->RsuuX[isq][iq][iNeut1]));
 514       fac += 4.0 * mHat * mf2 * sqrt(lambda) *
 515         real(coupSUSYPtr->LsuuX[isq][iq][iNeut1] * conj(coupSUSYPtr->RsuuX[isq][iq][iNeut1]));
 516     }
 517     fac *= 2.0/(1 - s2W);
 518   }
 519
 520   widNow = fac * preFac * ps * sqrt(lambda);
 521   if(DEBUG) {
 522     cout<<idRes<<":: id1:"<<id1Abs<<" id2:"<<id2Abs<<" Width: ";
 523     cout<<scientific<<widNow<<endl;
 524   }
 525   return;
 526 }
 527
 528 //==========================================================================
 529
 530 //  Class ResonanceChar
 531 //  Derived class for Neutralino Resonances
 532
 533 //--------------------------------------------------------------------------
 534
 535
 536 void ResonanceChar::initConstants(){
 537
 538   alpEM = coupSUSYPtr->alphaEM(mHat * mHat);
 539   s2W   = coupSUSYPtr->sin2W;
 540   return;
 541 }
 542
 543 //--------------------------------------------------------------------------
 544
 545 // Calculate various common prefactors for the current mass.
 546 void  ResonanceChar::calcPreFac(bool){
 547
 548   preFac = alpEM / (8.0 * s2W * pow(mHat,3));
 549   return;
 550 }
 551
 552 //--------------------------------------------------------------------------
 553
 554 // Calculate width for currently considered channel.
 555 void  ResonanceChar::calcWidth(bool){
 556
 557   widNow = 0.0;
 558   if(ps == 0.) return;
 559
 560   double fac = 0.0;
 561   kinFac = mHat * mHat - mf1 * mf1 + mf2 * mf2;
 562   kinFac2 = pow(mHat,4) + pow(mf1,4) - 2.0 * pow(mf2,4)
 563     + pow2(mHat) * pow2(mf2) + pow2(mf1) * pow2(mf2) - pow2(mHat) * pow2(mf1);
 564   lambda = lam(mHat*mHat , mf1*mf1 , mf2*mf2);
 565
 566   int idChar1 = typeChar(idRes);
 567   int idChar2 = typeChar(id1Abs);
 568   int idNeut1 = typeNeut(id1Abs);
 569
 570   if(idChar2>0 && id2Abs == 23){
 571     // ~chi_i -> chi_j + Z
 572     fac = kinFac2 * (norm(coupSUSYPtr->OLp[idChar1][idChar2])
 573                      + norm(coupSUSYPtr->ORp[idChar1][idChar2]));
 574     fac -= 12.0 * mHat * mf1 * pow2(mf2) *
 575       real(coupSUSYPtr->OLp[idChar1][idChar2]
 576            * conj(coupSUSYPtr->ORp[idChar1][idChar2]));
 577     fac /= pow2(mf2) * (1.0 - s2W);
 578   }
 579   else if(idNeut1>0 && id2Abs==24){
 580     // ~chi_i -> chi0_j + W- (or c.c.)
 581
 582     fac = kinFac2 * (norm(coupSUSYPtr->OL[idNeut1][idChar1]) + norm(coupSUSYPtr->OR[idNeut1][idChar1]));
 583     fac -= 12.0 * mHat * mf1 * pow2(mf2) *
 584       real(coupSUSYPtr->OL[idNeut1][idChar1] * conj(coupSUSYPtr->OR[idNeut1][idChar1]));
 585     fac /= pow2(mf2);
 586   }
 587   else if(id1Abs > 1000000 && id1Abs%100 < 7 && id2Abs < 7){
 588     // ~chi0_k -> ~q + q
 589     bool idown = (id1Abs%2 == 1);
 590     int iq = (id2Abs + 1 )/ 2;
 591     int isq = (abs(idRes)/1000000 == 2) ?
 592       (abs(idRes)%10+1)/2 + 3: (abs(idRes)%10+1)/2;
 593
 594     if(idown){
 595       fac = kinFac * (norm(coupSUSYPtr->LsduX[isq][iq][idChar1]) + norm(coupSUSYPtr->RsduX[isq][iq][idChar1]));
 596       fac += 4.0 * mHat * mf2 *
 597         real(coupSUSYPtr->LsduX[isq][iq][idChar1] * conj(coupSUSYPtr->RsduX[isq][iq][idChar1]));
 598     }
 599     else{
 600       fac = kinFac * (norm(coupSUSYPtr->LsudX[isq][iq][idChar1]) + norm(coupSUSYPtr->RsudX[isq][iq][idChar1]));
 601       fac += 4.0 * mHat * mf2 *
 602         real(coupSUSYPtr->LsudX[isq][iq][idChar1] * conj(coupSUSYPtr->RsudX[isq][iq][idChar1]));
 603     }
 604     fac *= 2.0/(1 - s2W);
 605   }
 606
 607   widNow = fac * preFac * ps * sqrt(lambda) ;
 608   if(DEBUG) {
 609     cout<<idRes<<":: id1:"<<id1Abs<<" id2:"<<id2Abs<<" Width: ";
 610     cout<<scientific<<widNow<<endl;
 611   }
 612   return;
 613 }
 614
 615
 616 //==========================================================================
 617
 618 //Return neutralino code; zero if not a neutralino
 619
 620 int SUSYResonanceWidths::typeNeut(int idPDG) {
 621   int type = 0;
 622   int idAbs = abs(idPDG);
 623   if(idAbs == 1000022) type = 1;
 624   else if(idAbs == 1000023) type = 2;
 625   else if(idAbs == 1000025) type = 3;
 626   else if(idAbs == 1000035) type = 4;
 627   else if(coupSUSYPtr->isNMSSM && idAbs == 1000045) type = 5;
 628   return type;
 629
 630 }
 631
 632
 633 //--------------------------------------------------------------------------
 634
 635 //Check whether particle is a Chargino
 636
 637 int SUSYResonanceWidths::typeChar(int idPDG) {
 638   int type = 0;
 639   if(abs(idPDG) == 1000024) type = 1;
 640   else if (abs(idPDG) == 1000037)type = 2;
 641   return type;
 642 }
 643
 644 //--------------------------------------------------------------------------
 645
 646 // Function for Kallen function
 647
 648 double SUSYResonanceWidths::lam(double x, double y, double z){
 649
 650   double val = x*x + y*y + z*z - 2.0* (x*y + y*z + z*x);
 651   return val;
 652
 653 }
 654
 655
 656 //--------------------------------------------------------------------------
 657
 658
 659
 660 } //end namespace Pythia8
 661
 662