PYTHIA8/pythia8170/src/SusyResonanceWidths.cxx

   1 // SusyResonanceWidths.cc is a part of the PYTHIA event generator.
   2 // Copyright (C) 2012 Torbjorn Sjostrand
   3 // Authors: 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 // SUSY Resonance widths classes.
   9
  10 #include "SusyResonanceWidths.h"
  11 #include "SusyCouplings.h"
  12 #include "PythiaComplex.h"
  13
  14 namespace Pythia8 {
  15
  16 //==========================================================================
  17
  18 // WidthFunctions Class
  19
  20 // Contains functions to be integrated for calculating the 3-body
  21 // decay widths
  22
  23 //--------------------------------------------------------------------------
  24
  25 void WidthFunction::init( ParticleData* particleDataPtrIn,
  26   CoupSUSY* coupSUSYPtrIn) {
  27
  28   particleDataPtr = particleDataPtrIn;
  29   coupSUSYPtr = coupSUSYPtrIn;
  30
  31 }
  32
  33 //--------------------------------------------------------------------------
  34
  35 void WidthFunction::setInternal2(int idResIn, int id1In, int id2In,
  36   int id3In, int idIntIn) {
  37
  38   // Res -> 1,2,3
  39   idRes = idResIn;
  40   id1 = id1In;
  41   id2 = id2In;
  42   id3 = id3In;
  43   idInt = idIntIn;
  44
  45   mRes = particleDataPtr->m0(idRes);
  46   m1 = particleDataPtr->m0(id1);
  47   m2 = particleDataPtr->m0(id2);
  48   m3 = particleDataPtr->m0(id3);
  49
  50   // Internal propagator
  51   mInt = particleDataPtr->m0(idInt);
  52   gammaInt = particleDataPtr->mWidth(idInt);
  53
  54   return;
  55 }
  56
  57 //--------------------------------------------------------------------------
  58
  59 double WidthFunction::function(double) {
  60
  61   cout<<"Warning using dummy width function"<<endl;
  62   return 0.;
  63 }
  64
  65 //--------------------------------------------------------------------------
  66
  67 double WidthFunction::function(double,double) {
  68
  69   cout<<"Warning using dummy width function"<<endl;
  70   return 0.;
  71 }
  72
  73 //==========================================================================
  74
  75 // Psi, Upsilon and Phi classes.
  76
  77 //--------------------------------------------------------------------------
  78
  79 void Psi::setInternal (int idResIn, int id1In, int id2In, int id3In,
  80   int idIntIn, int) {
  81
  82   setInternal2(idResIn, id1In, id2In, id3In, idIntIn);
  83
  84   mInt = particleDataPtr->m0(idInt);
  85   gammaInt = particleDataPtr->mWidth(idInt);
  86   iX = coupSUSYPtr->typeNeut(idRes);
  87   iQ = (id3+1)/2;
  88   iSq = (idInt>1000000)? 3 + (idInt%10+1)/2 :  (idInt%10+1)/2;
  89   isSqDown = (idInt % 2 == 1)? true : false;
  90   m1 = particleDataPtr->m0(id1);
  91   m2 = particleDataPtr->m0(id2);
  92   m3 = particleDataPtr->m0(id3);
  93   return;
  94 }
  95
  96 //--------------------------------------------------------------------------
  97
  98 void Upsilon::setInternal (int idResIn, int id1In, int id2In, int id3In,
  99   int idIntIn, int idInt2In) {
 100
 101   setInternal2(idResIn, id1In, id2In, id3In, idIntIn);
 102
 103   idInt2 = idInt2In;
 104   mInt = particleDataPtr->m0(idInt);
 105   gammaInt = particleDataPtr->mWidth(idInt);
 106   mInt2 = particleDataPtr->m0(idInt2);
 107   gammaInt2 = particleDataPtr->mWidth(idInt2);
 108
 109   iX = coupSUSYPtr->typeNeut(idRes);
 110   iQ = (id3+1)/2;
 111   iSq  = (idInt>1000000)? 3 + (idInt%10+1)/2 :  (idInt%10+1)/2;
 112   iSq2 = (idInt2>1000000)? 3 + (idInt2%10+1)/2 :  (idInt2%10+1)/2;
 113   isSqDown = (idIntIn % 2 == 1)? true : false;
 114   m1 = particleDataPtr->m0(id1);
 115   m2 = particleDataPtr->m0(id2);
 116   m3 = particleDataPtr->m0(id3);
 117   return;
 118 }
 119
 120 //--------------------------------------------------------------------------
 121
 122 void Phi::setInternal (int idResIn, int id1In, int id2In, int id3In,
 123   int idIntIn, int idInt2In) {
 124
 125   setInternal2(idResIn, id1In, id2In, id3In, idIntIn);
 126
 127   idInt2 = idInt2In;
 128   mInt = particleDataPtr->m0(idInt);
 129   gammaInt = particleDataPtr->mWidth(idInt);
 130   mInt2 = particleDataPtr->m0(idInt2);
 131   gammaInt2 = particleDataPtr->mWidth(idInt2);
 132
 133   iX = coupSUSYPtr->typeNeut(idRes);
 134   iQ = (id3+1)/2;
 135   iSq  = (idInt>1000000)? 3 + (idInt%10+1)/2 :  (idInt%10+1)/2;
 136   iSq2 = (idInt2>1000000)? 3 + (idInt2%10+1)/2 :  (idInt2%10+1)/2;
 137   isSqDown = (idIntIn % 2 == 1)? true : false;
 138   m1 = particleDataPtr->m0(id1);
 139   m2 = particleDataPtr->m0(id2);
 140   m3 = particleDataPtr->m0(id3);
 141   return;
 142 }
 143
 144 //--------------------------------------------------------------------------
 145
 146 double Psi::function(double m12sq) {
 147
 148   double R, factor1, factor2, value;
 149
 150   // Check that the propagators are offshell
 151   if(m12sq > pow2(mInt) || abs(m12sq - pow2(mInt)) < gammaInt) return 0;
 152
 153   R = 1.0/(pow2(m12sq-pow2(mInt)) + pow2(gammaInt*mInt));
 154   if(isSqDown){
 155     factor1 = (norm(coupSUSYPtr->LsddX[iSq][iQ][iX])
 156                + norm(coupSUSYPtr->RsddX[iSq][iQ][iX]))*
 157       (mRes*mRes + m3*m3 - m12sq);
 158     factor2 = 4.0 * real(coupSUSYPtr->LsddX[iSq][iQ][iX]
 159                          * conj(coupSUSYPtr->RsddX[iSq][iQ][iX])) * m3 * mRes;
 160
 161   } else {
 162     factor1 = (norm(coupSUSYPtr->LsuuX[iSq][iQ][iX])
 163                + norm(coupSUSYPtr->RsuuX[iSq][iQ][iX]))*
 164       (mRes*mRes + m3*m3 - m12sq);
 165     factor2 = 4.0 * real(coupSUSYPtr->LsuuX[iSq][iQ][iX]
 166                          * conj(coupSUSYPtr->RsuuX[iSq][iQ][iX])) * m3 * mRes;
 167   }
 168
 169   value = R * (m12sq - m1*m1 - m2*m2) * (factor1+factor2);
 170
 171   return value;
 172 }
 173
 174
 175 //--------------------------------------------------------------------------
 176
 177 double Upsilon::function(double m12sq) {
 178
 179   double R1,R2, S, factor1, factor2, value;
 180
 181   // Check that the propagators are offshell
 182   if(m12sq > pow2(mInt) || abs(m12sq - pow2(mInt)) < gammaInt) return 0;
 183   if(m12sq > pow2(mInt2) || abs(m12sq - pow2(mInt2)) < gammaInt2) return 0;
 184
 185   R1 = 1.0/(pow2(m12sq-pow2(mInt)) + pow2(gammaInt*mInt));
 186   R2 =  1.0/(pow2(m12sq-pow2(mInt2)) + pow2(gammaInt2*mInt2));
 187   S = R1 * R2 * ((m12sq - pow2(mInt)) * (m12sq - pow2(mInt2))
 188                  + gammaInt * mInt * gammaInt2 * mInt2);
 189
 190   if(isSqDown){
 191     factor1 = real(coupSUSYPtr->LsddX[iSq][iQ][iX]
 192                    * conj(coupSUSYPtr->LsddX[iSq2][iQ][iX]))
 193       + real(coupSUSYPtr->RsddX[iSq][iQ][iX]
 194              * conj(coupSUSYPtr->RsddX[iSq2][iQ][iX]));
 195     factor2 = real(coupSUSYPtr->LsddX[iSq][iQ][iX]
 196                    * conj(coupSUSYPtr->RsddX[iSq2][iQ][iX]))
 197       + real(coupSUSYPtr->RsddX[iSq][iQ][iX]
 198              * conj(coupSUSYPtr->LsddX[iSq2][iQ][iX]));
 199   }else{
 200     factor1 = real(coupSUSYPtr->LsuuX[iSq][iQ][iX]
 201                    * conj(coupSUSYPtr->LsuuX[iSq2][iQ][iX]))
 202       + real(coupSUSYPtr->RsuuX[iSq][iQ][iX]
 203              * conj(coupSUSYPtr->RsuuX[iSq2][iQ][iX]));
 204     factor2 = real(coupSUSYPtr->LsuuX[iSq][iQ][iX]
 205                    * conj(coupSUSYPtr->RsuuX[iSq2][iQ][iX]))
 206       + real(coupSUSYPtr->RsuuX[iSq][iQ][iX]
 207              * conj(coupSUSYPtr->LsuuX[iSq2][iQ][iX]));
 208   }
 209
 210   value = S * (m12sq - pow2(m1) - pow2(m2)) *
 211     ( factor1 * (pow2(mRes) + pow2(m3) - m12sq) + 2.0 * factor2 * m3 * mRes);
 212
 213   //  cout<<"I1: "<<idInt<<" I2:"<<idInt2<<" factor1: "<<factor1
 214   //        <<" factor2:"<<factor2<<" value:"<<value<<endl;
 215
 216   return value;
 217
 218 }
 219
 220 //--------------------------------------------------------------------------
 221
 222 double Phi::function(double m12sqIn) {
 223
 224   m12sq =  m12sqIn;
 225   // Check that the propagators are offshell
 226   if(m12sq > pow2(mInt) || abs(m12sq - pow2(mInt)) < gammaInt) return 0;
 227
 228   double m23max, m23min, E2, E3;
 229
 230   E2 = (m12sq - pow2(m1) - pow2(m2))/(2.0 * sqrt(m12sq));
 231   E3 = (pow2(mRes) - m12sq - pow2(m3))/(2.0 * sqrt(m12sq));
 232   m23max = pow2(E2+E3) - (sqrt(E2*E2 - m2*m2) - sqrt(E3*E3 - m3*m3)) ;
 233   m23min = pow2(E2+E3) - (sqrt(E2*E2 - m2*m2) + sqrt(E3*E3 - m3*m3)) ;
 234
 235   if(E2 < m2 || E3 < m3){
 236     cout <<"Error in Phi:function"<<endl;
 237     return 0.;
 238   }
 239
 240   double integral2 = integrateGauss(m23min,m23max,1.0e-4);
 241   return integral2;
 242 }
 243
 244 //--------------------------------------------------------------------------
 245
 246 double Phi::function2(double m23sq) {
 247
 248   // Check that the propagators are offshell
 249   if(m23sq > pow2(mInt2) || abs(m23sq - pow2(mInt2)) < gammaInt2) return 0;
 250
 251   double R1, R2, S, factor1, factor2, factor3, factor4, value, fac;
 252   int iQ2;
 253
 254   R1 = 1.0/(pow2(m12sq-pow2(mInt)) + pow2(gammaInt*mInt));
 255   R2 = 1.0/(pow2(m12sq-pow2(mInt2)) + pow2(gammaInt2*mInt2));
 256   S = R1 * R2 * ((m12sq - pow2(mInt)) * (m12sq - pow2(mInt2))
 257     + gammaInt * mInt * gammaInt2 * mInt2);
 258
 259   fac = 1.0;
 260
 261   if(isSqDown){
 262     // Only factor is when both d_i and d_j are near.
 263     iQ2 = (id1%2 == 1)? (id1+1)/2 : (id2+1)/2;
 264
 265     if(mRes > mInt2 + particleDataPtr->m0(iQ2)) fac = 0.;
 266
 267     factor1 = m1 * m3 * real(coupSUSYPtr->LsddX[iSq][iQ][iX]
 268                              * conj(coupSUSYPtr->LsddX[iSq2][iQ2][iX]))
 269       * (m12sq + m23sq - pow2(m1) - pow2(m3));
 270     factor2 = m1 * mRes * real(coupSUSYPtr->RsddX[iSq][iQ][iX]
 271                                * conj(coupSUSYPtr->LsddX[iSq2][iQ2][iX]))
 272       * (m23sq - pow2(m2) - pow2(m3));
 273     factor3 = m3 * mRes * real(coupSUSYPtr->LsddX[iSq][iQ][iX]
 274                                * conj(coupSUSYPtr->RsddX[iSq2][iQ2][iX]))
 275       * (m12sq - pow2(m1) - pow2(m2));
 276     factor4 = m3 * mRes * real(coupSUSYPtr->RsddX[iSq][iQ][iX]
 277                                * conj(coupSUSYPtr->RsddX[iSq2][iQ2][iX]))
 278       * (m12sq * m23sq - pow2(m1 * m3) - pow2(m2 * mRes));
 279
 280   }else{
 281     // Factor A: u and d_1
 282     iQ2 = (id1+1)/2;
 283
 284     if(mRes > mInt2 + particleDataPtr->m0(iQ2)) fac = 0.;
 285
 286     factor1 = m1 * m3 * real(coupSUSYPtr->LsuuX[iSq][iQ][iX]
 287                              * conj(coupSUSYPtr->LsddX[iSq2][iQ2][iX]))
 288       * (m12sq + m23sq - pow2(m1) - pow2(m3));
 289     factor2 = m1 * mRes * real(coupSUSYPtr->RsuuX[iSq][iQ][iX]
 290                                * conj(coupSUSYPtr->LsddX[iSq2][iQ2][iX]))
 291       * (m23sq - pow2(m2) - pow2(m3));
 292     factor3 = m3 * mRes * real(coupSUSYPtr->LsuuX[iSq][iQ][iX]
 293                                * conj(coupSUSYPtr->RsddX[iSq2][iQ2][iX]))
 294       * (m12sq - pow2(m1) - pow2(m2));
 295     factor4 = m3 * mRes * real(coupSUSYPtr->RsuuX[iSq][iQ][iX]
 296                                * conj(coupSUSYPtr->RsddX[iSq2][iQ2][iX]))
 297       * (m12sq * m23sq - pow2(m1 * m3) - pow2(m2 * mRes));
 298
 299
 300     // Factor B: u and d_2; change 1 <=> 2
 301     iQ2 = (id2+1)/2;
 302
 303     if(mRes > mInt2 + particleDataPtr->m0(iQ2)) fac = 0.;
 304
 305     factor1 += m2 * m3 * real(coupSUSYPtr->LsuuX[iSq][iQ][iX]
 306                               * conj(coupSUSYPtr->LsddX[iSq2][iQ2][iX]))
 307       * (m12sq + m23sq - pow2(m2) - pow2(m3));
 308     factor2 += m2 * mRes * real(coupSUSYPtr->RsuuX[iSq][iQ][iX]
 309                                 * conj(coupSUSYPtr->LsddX[iSq2][iQ2][iX]))
 310       * (m23sq - pow2(m1) - pow2(m3));
 311     factor3 += m3 * mRes * real(coupSUSYPtr->LsuuX[iSq][iQ][iX]
 312                                 * conj(coupSUSYPtr->RsddX[iSq2][iQ2][iX]))
 313       * (m12sq - pow2(m2) - pow2(m1));
 314     factor4 += m3 * mRes * real(coupSUSYPtr->RsuuX[iSq][iQ][iX]
 315                                 * conj(coupSUSYPtr->RsddX[iSq2][iQ2][iX]))
 316       * (m12sq * m23sq - pow2(m2 * m3) - pow2(m1 * mRes));
 317   }
 318
 319   value = S * (factor1 + factor2 + factor3 + factor4);
 320
 321   //  cout<<"I1: "<<idInt<<" I2:"<<idInt2<<" factor1: "<<factor1
 322   //      <<" factor2:"<<factor2<<" value:"<<value<<endl;
 323
 324   return (fac * value);
 325 }
 326
 327 //--------------------------------------------------------------------------
 328
 329 double Phi::integrateGauss(double xlo, double xhi, double tol) {
 330
 331   //8-point unweighted
 332   static double x8[4]={0.96028985649753623,
 333                        0.79666647741362674,
 334                        0.52553240991632899,
 335                        0.18343464249564980};
 336   static double w8[4]={0.10122853629037626,
 337                        0.22238103445337447,
 338                        0.31370664587788729,
 339                        0.36268378337836198};
 340   //16-point unweighted
 341   static double x16[8]={0.98940093499164993,
 342                        0.94457502307323258,
 343                        0.86563120238783174,
 344                        0.75540440835500303,
 345                        0.61787624440264375,
 346                        0.45801677765722739,
 347                        0.28160355077925891,
 348                        0.09501250983763744};
 349   static double w16[8]={0.027152459411754095,
 350                        0.062253523938647893,
 351                        0.095158511682492785,
 352                        0.12462897125553387,
 353                        0.14959598881657673,
 354                        0.16915651939500254,
 355                        0.18260341504492359,
 356                        0.18945061045506850};
 357   //boundary checks
 358   if (xlo == xhi) {
 359     cerr<<"xlo = xhi"<<endl;
 360     return 0.0;
 361   }
 362   if ( xlo > xhi ) {
 363     cerr<<" (integrateGauss:) -> xhi < xlo"<<endl;
 364     return 0.0;
 365   }
 366   //initialize
 367   double sum = 0.0;
 368   double c = 0.001/abs(xhi-xlo);
 369   double zlo = xlo;
 370   double zhi = xhi;
 371
 372   bool nextbin = true;
 373
 374   while ( nextbin ) {
 375
 376     double zmi = 0.5*(zhi+zlo); // midpoint
 377     double zmr = 0.5*(zhi-zlo); // midpoint, relative to zlo
 378
 379     //calculate 8-point and 16-point quadratures
 380     double s8=0.0;
 381     for (int i=0;i<4;i++) {
 382       double dz = zmr * x8[i];
 383       s8 += w8[i]*(function2(zmi+dz) + function2(zmi-dz));
 384     }
 385     s8 *= zmr;
 386     double s16=0.0;
 387     for (int i=0;i<8;i++) {
 388       double dz = zmr * x16[i];
 389       s16 += w16[i]*(function2(zmi+dz) + function2(zmi-dz));
 390     }
 391     s16 *= zmr;
 392     if (abs(s16-s8) < tol*(1+abs(s16))) {
 393       //precision in this bin OK, add to cumulative and go to next
 394       nextbin=true;
 395       sum += s16;
 396       //next bin: LO = end of current, HI = end of integration region.
 397       zlo=zhi;
 398       zhi=xhi;
 399       if ( zlo == zhi ) nextbin = false;
 400     } else {
 401       //precision in this bin not OK, subdivide.
 402       if (1.0 + c*abs(zmr) == 1.0) {
 403         cerr << " (integrateGauss:) too high accuracy required"<<endl;
 404         sum = 0.0 ;
 405         break;
 406       }
 407       zhi=zmi;
 408       nextbin=true;
 409     }
 410   }
 411   return sum;
 412 }
 413
 414 //==========================================================================
 415
 416 // The SUSYResonanceWidths Class
 417 // Derived class for SUSY resonances
 418
 419 const bool SUSYResonanceWidths::DEBUG = false;
 420
 421 //--------------------------------------------------------------------------
 422
 423 bool SUSYResonanceWidths::init(Info* infoPtrIn, Settings* settingsPtrIn,
 424    ParticleData* particleDataPtrIn, Couplings* couplingsPtrIn) {
 425
 426   // Save pointers.
 427   infoPtr         = infoPtrIn;
 428   settingsPtr     = settingsPtrIn;
 429   particleDataPtr = particleDataPtrIn;
 430   coupSUSYPtr     = (couplingsPtrIn->isSUSY ? (CoupSUSY *) couplingsPtrIn: 0 );
 431
 432   // No width calculation necessary for SM-only
 433   bool calcWidthAllow = true;
 434   if(!couplingsPtrIn->isSUSY) calcWidthAllow = false;
 435
 436   // Minimal decaying-resonance width. Minimal phase space for meMode = 103.
 437   minWidth     = settingsPtr->parm("ResonanceWidths:minWidth");
 438   minThreshold = settingsPtr->parm("ResonanceWidths:minThreshold");
 439
 440   // Pointer to particle species.
 441   particlePtr  = particleDataPtr->particleDataEntryPtr(idRes);
 442   if (particlePtr == 0) {
 443     infoPtr->errorMsg("Error in ResonanceWidths::init:"
 444       " unknown resonance identity code");
 445     return false;
 446   }
 447
 448   // Generic particles should not have meMode < 100, but allow
 449   // some exceptions: not used Higgses and not used Technicolor.
 450   if (idRes == 35 || idRes == 36 || idRes == 37
 451     || idRes/1000000 == 3) isGeneric = false;
 452
 453   // Resonance properties: antiparticle, mass, width
 454   hasAntiRes   = particlePtr->hasAnti();
 455   mRes         = particlePtr->m0();
 456   GammaRes     = particlePtr->mWidth();
 457   m2Res        = mRes*mRes;
 458
 459   // For very narrow resonances assign fictitious small width.
 460   if (GammaRes < minWidth) GammaRes = 0.1 * minWidth;
 461   GamMRat      = GammaRes / mRes;
 462
 463   // Secondary decay chains by default all on.
 464   openPos      = 1.;
 465   openNeg      = 1.;
 466
 467   // Allow option where on-shell width is forced to current value.
 468   doForceWidth = particlePtr->doForceWidth();
 469   forceFactor  = 1.;
 470
 471   // Check that resonance OK.
 472   if (particlePtr == 0) infoPtr->errorMsg("Error in ResonanceWidths::init:"
 473       " unknown resonance identity code");
 474
 475   // Check if decay table was read in via SLHA
 476   bool hasDecayTable = false;
 477   if(calcWidthAllow) {
 478     for(unsigned int iDec = 1; iDec < (coupSUSYPtr->slhaPtr)->decays.size();
 479       iDec++)
 480       if(!hasDecayTable) hasDecayTable
 481         = ((coupSUSYPtr->slhaPtr)->decays[iDec].getId() == abs(idRes));
 482   }
 483
 484   if(hasDecayTable && settingsPtr->flag("SLHA:useDecayTable")) {
 485     calcWidthAllow = false;
 486     if(DEBUG) cout<<"Found decay table for:"<<idRes<<endl;
 487   }
 488
 489   // Initialize constants used for a resonance.
 490   if(calcWidthAllow) {
 491     mHat          = mRes;
 492     initConstants();
 493     calcPreFac(true);
 494   }
 495
 496   // Reset quantities to sum. Declare variables inside loop.
 497   double widTot = 0.;
 498   double widPos = 0.;
 499   double widNeg = 0.;
 500   int    idNow, idAnti;
 501   double openSecPos, openSecNeg;
 502
 503   // Loop over all decay channels. Basic properties of channel.
 504   for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
 505     iChannel    = i;
 506     onMode      = particlePtr->channel(i).onMode();
 507     meMode      = particlePtr->channel(i).meMode();
 508     mult        = particlePtr->channel(i).multiplicity();
 509     widNow      = 0.;
 510
 511     // Warn if not relevant meMode.
 512     if ( meMode < 0 || meMode > 103 || (isGeneric && meMode < 100) ) {
 513       infoPtr->errorMsg("Error in ResonanceWidths::init:"
 514         " resonance meMode not acceptable");
 515     }
 516
 517     // Calculation of SUSY particle widths
 518     if (meMode == 103 && GammaRes > 0. && calcWidthAllow &&
 519         (!settingsPtr->flag("SLHA:useDecayTable") || !hasDecayTable)) {
 520       // Read out information on channel: primarily use first two.
 521       id1       = particlePtr->channel(i).product(0);
 522       id2       = particlePtr->channel(i).product(1);
 523       id1Abs    = abs(id1);
 524       id2Abs    = abs(id2);
 525
 526       // Order first two in descending order of absolute values.
 527       if (id2Abs > id1Abs) {swap( id1, id2); swap( id1Abs, id2Abs);}
 528
 529       // Allow for third product to be treated in derived classes.
 530       if (mult > 2) {
 531         id3     = particlePtr->channel(i).product(2);
 532         id3Abs  = abs(id3);
 533
 534         // Also order third into descending order of absolute values.
 535         if (id3Abs > id2Abs) {swap( id2, id3); swap( id2Abs, id3Abs);}
 536         if (id2Abs > id1Abs) {swap( id1, id2); swap( id1Abs, id2Abs);}
 537       }
 538
 539       // Read out masses. Calculate two-body phase space.
 540       mf1       = particleDataPtr->m0(id1Abs);
 541       mf2       = particleDataPtr->m0(id2Abs);
 542       mr1       = pow2(mf1 / mHat);
 543       mr2       = pow2(mf2 / mHat);
 544       ps        = (mHat < mf1 + mf2 + MASSMARGIN) ? 0.
 545         : sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 ) * pow2(mHat);
 546       if (mult > 2) {
 547         mf3     = particleDataPtr->m0(id3Abs);
 548         mr3     = pow2(mf3 / mHat);
 549
 550         //Check phase space
 551         ps = 1.0;
 552         if(mHat < mf1 + mf2 + mf3 + MASSMARGIN) ps = 0.;
 553       }
 554
 555       // Let derived class calculate width for channel provided.
 556       calcWidth(true);
 557     }
 558
 559     // Width calculated based on stored values.
 560     else widNow = GammaRes * particlePtr->channel(i).bRatio();
 561
 562     // Find secondary open fractions of partial width.
 563     openSecPos  = 1.;
 564     openSecNeg  = 1.;
 565     if (widNow > 0.) for (int j = 0; j < mult; ++j) {
 566       idNow     = particlePtr->channel(i).product(j);
 567       idAnti    = (particleDataPtr->hasAnti(idNow)) ? -idNow : idNow;
 568       // Secondary widths not yet initialized for heavier states,
 569       // so have to assume unit open fraction there.
 570       if (idNow == 23 || abs(idNow) == 24
 571         || particleDataPtr->m0(abs(idNow)) < mRes) {
 572         openSecPos *= particleDataPtr->resOpenFrac(idNow);
 573         openSecNeg *= particleDataPtr->resOpenFrac(idAnti);
 574       }
 575     }
 576
 577     // Store partial widths and secondary open fractions.
 578     particlePtr->channel(i).onShellWidth(widNow);
 579     particlePtr->channel(i).openSec( idRes, openSecPos);
 580     particlePtr->channel(i).openSec(-idRes, openSecNeg);
 581
 582     // Update sum over all channnels and over open channels only.
 583     widTot     += widNow;
 584     if (onMode == 1 || onMode == 2) widPos += widNow * openSecPos;
 585     if (onMode == 1 || onMode == 3) widNeg += widNow * openSecNeg;
 586   }
 587
 588   // If no decay channels are open then set particle stable and done.
 589   if (widTot < minWidth) {
 590     particlePtr->setMayDecay(false, false);
 591     particlePtr->setMWidth(0., false);
 592     for (int i = 0; i < particlePtr->sizeChannels(); ++i)
 593       particlePtr->channel(i).bRatio( 0., false);
 594     return true;
 595   }
 596
 597   // Normalize branching ratios to unity.
 598   double bRatio;
 599   for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
 600     bRatio      = particlePtr->channel(i).onShellWidth() / widTot;
 601     particlePtr->channel(i).bRatio( bRatio, false);
 602   }
 603
 604   // Optionally force total width by rescaling of all partial ones.
 605   if (doForceWidth) {
 606     forceFactor = GammaRes / widTot;
 607     for (int i = 0; i < particlePtr->sizeChannels(); ++i)
 608       particlePtr->channel(i).onShellWidthFactor( forceFactor);
 609   }
 610
 611   // Else update newly calculated partial width.
 612   else {
 613     particlePtr->setMWidth(widTot, false);
 614     GammaRes    = widTot;
 615   }
 616
 617   // Updated width-to-mass ratio. Secondary widths for open.
 618   GamMRat       = GammaRes / mRes;
 619   openPos       = widPos / widTot;
 620   openNeg       = widNeg / widTot;
 621
 622   // Done.
 623   return true;
 624
 625 }
 626
 627 //==========================================================================
 628
 629 // The ResonanceSquark class
 630 // Derived class for Squark resonances
 631
 632 //--------------------------------------------------------------------------
 633
 634 // Initialize constants.
 635
 636 void ResonanceSquark::initConstants() {
 637
 638   // Locally stored properties and couplings.
 639   alpS  = coupSUSYPtr->alphaS(mHat * mHat );
 640   alpEM = coupSUSYPtr->alphaEM(mHat * mHat);
 641   s2W   = coupSUSYPtr->sin2W;
 642 }
 643
 644 //--------------------------------------------------------------------------
 645
 646 // Calculate various common prefactors for the current mass.
 647
 648 void ResonanceSquark::calcPreFac(bool) {
 649
 650   // Common coupling factors.
 651   preFac = 1.0 / (s2W * pow(mHat,3));
 652
 653 }
 654
 655 //--------------------------------------------------------------------------
 656
 657 // Calculate width for currently considered channel.
 658
 659 void ResonanceSquark::calcWidth(bool) {
 660
 661   // Squark type -- in u_i/d_i and generation
 662   int ksusy = 1000000;
 663   bool idown = (abs(idRes)%2 == 0 ? false : true);
 664   int isq = (abs(idRes)/ksusy == 2) ?
 665     (abs(idRes)%10+1)/2 + 3: (abs(idRes)%10+1)/2;
 666   // int isqgen = (abs(idRes)%10 + 1)/2;
 667
 668   // Check that mass is above threshold.
 669   if(ps == 0.) return;
 670   else{
 671     // Two-body decays
 672     kinFac = (mHat * mHat - mf1 * mf1 - mf2 * mf2);
 673
 674     double fac = 0.0 , wid = 0.0;
 675
 676     //RPV decays
 677     //Case 1a:  UDD-type
 678     if(id1Abs < 7 && id2Abs < 7){
 679
 680       // Quark generations
 681       int iq1 = (id1Abs + 1)/2;
 682       int iq2 = (id2Abs + 1)/2;
 683
 684       // Check for RPV UDD couplings
 685       if(!coupSUSYPtr->isUDD) {widNow = 0; return;}
 686
 687       // ~q -> q_i + q_j
 688
 689       fac = 2.0 * kinFac / (16.0 * M_PI * pow(mHat,3));
 690       wid = 0.0;
 691       if(idown) {
 692         if ((id1Abs+id2Abs)%2 == 1){
 693           if(id1Abs%2==1)
 694             for(int isq2 = 1; isq2 < 4; isq2++)
 695               wid += norm(coupSUSYPtr->rvUDD[iq2][iq1][isq2]
 696                    * coupSUSYPtr->Rdsq[isq][isq2+3]);
 697           else
 698             for(int isq2 = 1; isq2 < 4; isq2++)
 699               wid += norm(coupSUSYPtr->rvUDD[iq1][iq2][isq2]
 700                    * coupSUSYPtr->Rdsq[isq][isq2+3]);
 701         }
 702       }
 703       else {
 704         if ((id1Abs+id2Abs)%2 != 0) widNow = 0.0;
 705         else
 706           for(int isq2 = 1; isq2 < 4; isq2++)
 707             wid += norm(coupSUSYPtr->rvUDD[isq2][iq1][iq2]
 708                  * coupSUSYPtr->Rusq[isq][isq2+3]);
 709       }
 710   }
 711
 712     //Case 1b:  LQD-type
 713     else if(id1Abs < 17 && id2Abs < 7){
 714       if(!coupSUSYPtr->isLQD) {widNow = 0; return;}
 715
 716       int ilep = (id1Abs - 9)/2;
 717       int iq = (id2Abs + 1)/2;
 718
 719       fac = kinFac / (16.0 * M_PI * pow(mHat,3));
 720       wid = 0.0;
 721       if(idown){
 722         if(iq%2 == 0){
 723           // q is up-type; ~q is right-handed down type
 724           for(int isq2=1; isq2<3; isq2++)
 725             wid += norm(coupSUSYPtr->Rdsq[isq][isq2+3]
 726                  * coupSUSYPtr->rvLQD[ilep][iq][isq2]);
 727         }else{
 728           //q is down type; ~q left-handed down-type
 729           for(int isq2=1; isq2<3; isq2++)
 730             wid += norm(coupSUSYPtr->Rdsq[isq][isq2]
 731                  * coupSUSYPtr->rvLQD[ilep][isq2][isq2]);
 732         }
 733       }
 734       else{
 735         if(iq%2 == 0) {widNow = 0.0; return;}
 736         // q is down type; ~q is left-handed up-type
 737         for(int isq2=1; isq2<3; isq2++)
 738           wid += norm(coupSUSYPtr->Rusq[isq][isq2]
 739                * coupSUSYPtr->rvLQD[ilep][isq2][iq]);
 740       }
 741     }
 742
 743     //Case 2: quark + gaugino
 744     else if (id1Abs > ksusy && id2Abs < 7) {
 745
 746       int iq = (id2Abs + 1)/2;
 747
 748       // ~q -> ~g + q
 749       if(id1Abs == 1000021 && idRes%10 == id2Abs) {
 750         // Removed factor of s2W in denominator: strong process -- no EW
 751         fac = 2.0 * alpS / (3.0 * pow3(mHat));
 752         if(idown)
 753           wid = kinFac * (norm(coupSUSYPtr->LsddG[isq][iq])
 754               + norm(coupSUSYPtr->RsddG[isq][iq]))
 755               - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsddG[isq][iq]
 756               * conj(coupSUSYPtr->RsddG[isq][iq]));
 757         else
 758           wid = kinFac * (norm(coupSUSYPtr->LsuuG[isq][iq])
 759               + norm(coupSUSYPtr->RsuuG[isq][iq]))
 760               - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsuuG[isq][iq]
 761               * conj(coupSUSYPtr->RsuuG[isq][iq]));
 762       }
 763       else
 764         for(int i=1; i<6 ; i++){
 765           // ~q -> ~chi0 + q
 766           if(coupSUSYPtr->idNeut(i)==id1Abs && idRes%2 == id2Abs%2){
 767             fac = alpEM *  preFac / (2.0 * (1 - s2W));
 768             if(idown)
 769               wid = kinFac * (norm(coupSUSYPtr->LsddX[isq][iq][i])
 770                   + norm(coupSUSYPtr->RsddX[isq][iq][i]))
 771                   - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsddX[isq][iq][i]
 772                   * conj(coupSUSYPtr->RsddX[isq][iq][i]));
 773             else
 774               wid = kinFac * (norm(coupSUSYPtr->LsuuX[isq][iq][i])
 775                   + norm(coupSUSYPtr->RsuuX[isq][iq][i]))
 776                   - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsuuX[isq][iq][i]
 777                   * conj(coupSUSYPtr->RsuuX[isq][iq][i]));
 778           }
 779
 780           // ~q -> chi- + q
 781           else if (i < 3 && coupSUSYPtr->idChar(i)==id1Abs
 782             && idRes%2 != id2Abs%2){
 783
 784             fac = alpEM *  preFac / (4.0 * (1 - s2W));
 785             if(idown)
 786               wid = kinFac * (norm(coupSUSYPtr->LsduX[isq][iq][i])
 787                   + norm(coupSUSYPtr->RsduX[isq][iq][i]))
 788                   - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsduX[isq][iq][i]
 789                   * conj(coupSUSYPtr->RsduX[isq][iq][i]));
 790             else
 791               wid = kinFac * (norm(coupSUSYPtr->LsudX[isq][iq][i])
 792                   + norm(coupSUSYPtr->RsudX[isq][iq][i]))
 793                   - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsudX[isq][iq][i]
 794                   * conj(coupSUSYPtr->RsudX[isq][iq][i]));
 795           }
 796         }
 797     }
 798
 799     //Case 3: ~q_i -> ~q_j + Z/W
 800     else if (id1Abs > ksusy && id1Abs%100 < 7
 801       && (id2Abs == 23 || id2Abs == 24)){
 802
 803       // factor of lambda^(3/2) = ps^(3/2) ;
 804       fac = alpEM * preFac/(16.0 * pow2(particleDataPtr->m0(id2Abs))
 805           * (1.0 - s2W)) * pow2(ps) ;
 806
 807       int isq2 = (id1Abs/ksusy == 2) ? (id1Abs%10+1)/2 + 3: (id1Abs%10+1)/2;
 808
 809       if(id2Abs == 23 && id1Abs%2 == idRes%2){
 810         if(idown)
 811           wid = norm(coupSUSYPtr->LsdsdZ[isq][isq2]
 812                      + coupSUSYPtr->RsdsdZ[isq][isq2]);
 813         else
 814           wid = norm(coupSUSYPtr->LsusuZ[isq][isq2]
 815                      + coupSUSYPtr->RsusuZ[isq][isq2]);
 816       }
 817       else if (id2Abs == 24 && id1Abs%2 != idRes%2){
 818         if(idown)
 819           wid = norm(coupSUSYPtr->LsusdW[isq2][isq]);
 820         else
 821           wid = norm(coupSUSYPtr->LsusdW[isq][isq2]);
 822       }
 823     }
 824
 825     // TODO: Case ~q_i -> ~q_j + h/H
 826     widNow = fac * wid * ps;
 827     if(DEBUG) cout<<idRes<<":: id1:"<<id1Abs<<" id2:"<<id2Abs
 828                   <<" Width: "<<widNow<<endl;
 829     return;
 830   }
 831
 832 }
 833
 834 //==========================================================================
 835
 836 // The ResonanceGluino class
 837 // Derived class for Gluino resonances
 838
 839 //--------------------------------------------------------------------------
 840
 841 // Initialize constants.
 842
 843 void ResonanceGluino::initConstants() {
 844
 845   // Locally stored properties and couplings.
 846   alpS  = coupSUSYPtr->alphaS(mHat * mHat);
 847   return;
 848 }
 849
 850 //--------------------------------------------------------------------------
 851
 852 // Calculate various common prefactors for the current mass.
 853
 854 void ResonanceGluino::calcPreFac(bool) {
 855   // Common coupling factors.
 856   preFac = alpS /( 8.0 * pow(mHat,3));
 857   return;
 858 }
 859
 860 //--------------------------------------------------------------------------
 861
 862 // Calculate width for currently considered channel.
 863
 864 void ResonanceGluino::calcWidth(bool) {
 865
 866
 867   widNow = 0.0;
 868   if(ps == 0.) return;
 869   kinFac = (mHat * mHat - mf1 * mf1 + mf2 * mf2);
 870
 871   if(id1Abs > 1000000 && (id1Abs % 100) < 7 && id2Abs < 7) {
 872
 873     int isq = (abs(id1Abs)/1000000 == 2) ? (abs(id1Abs)%10+1)/2 + 3
 874                                          : (abs(id1Abs)%10+1)/2;
 875     bool idown = id2Abs%2;
 876     int iq = (id2Abs + 1)/2;
 877
 878     // ~g -> ~q + q
 879     if(idown){
 880       widNow = kinFac * (norm(coupSUSYPtr->LsddG[isq][iq])
 881              + norm(coupSUSYPtr->RsddG[isq][iq]))
 882              + 4.0 * mHat * mf2 * real(coupSUSYPtr->LsddG[isq][iq]
 883              * conj(coupSUSYPtr->RsddG[isq][iq]));
 884     }
 885     else{
 886       widNow = kinFac * (norm(coupSUSYPtr->LsuuG[isq][iq])
 887              + norm(coupSUSYPtr->RsuuG[isq][iq]))
 888              + 4.0 * mHat * mf2 * real(coupSUSYPtr->LsuuG[isq][iq]
 889              * conj(coupSUSYPtr->RsuuG[isq][iq]));
 890     }
 891     widNow = widNow * preFac * ps;
 892     if(DEBUG) {
 893       cout<<"Gluino:: id1:"<<id1Abs<<" id2:"<<id2Abs<<" Width: ";
 894       cout<<scientific<<widNow<<endl;
 895     }
 896     return;
 897   }
 898 }
 899
 900 //==========================================================================
 901
 902 //  Class ResonanceNeut
 903 //  Derived class for Neutralino Resonances
 904 //
 905 //--------------------------------------------------------------------------
 906
 907
 908 void ResonanceNeut::initConstants(){
 909
 910   alpEM = coupSUSYPtr->alphaEM(mHat * mHat);
 911   s2W   = coupSUSYPtr->sin2W;
 912
 913   // Initialize functions for calculating 3-body widths
 914   psi.init(particleDataPtr,coupSUSYPtr);
 915   phi.init(particleDataPtr,coupSUSYPtr);
 916   upsil.init(particleDataPtr,coupSUSYPtr);
 917 }
 918
 919 //--------------------------------------------------------------------------
 920
 921 // Calculate various common prefactors for the current mass.
 922 void  ResonanceNeut::calcPreFac(bool){
 923
 924   // Common coupling factors.
 925   preFac = alpEM / (8.0 * s2W * pow(mHat,3));
 926   return;
 927 }
 928
 929 //--------------------------------------------------------------------------
 930
 931 // Calculate width for currently considered channel.
 932 void  ResonanceNeut::calcWidth(bool){
 933
 934   widNow = 0.0;
 935
 936   if(ps ==0.) return;
 937   else if(mult ==2){
 938     // Two-body decays
 939
 940     kinFac = mHat * mHat - mf1 * mf1 + mf2 * mf2;
 941     kinFac2 = pow(mHat,4) + pow(mf1,4) - 2.0 * pow(mf2,4)
 942             + pow2(mHat) * pow2(mf2) + pow2(mf1) * pow2(mf2)
 943             - 2.0 * pow2(mHat) * pow2(mf1);
 944
 945     // Stable lightest neutralino
 946     if(idRes == 1000022) return;
 947
 948     double fac = 0.0;
 949     int iNeut1 = coupSUSYPtr->typeNeut(idRes);
 950     int iNeut2 = coupSUSYPtr->typeNeut(id1Abs);
 951     int iChar1 = coupSUSYPtr->typeChar(id1Abs);
 952
 953     if(iNeut2>0 && id2Abs == 23){
 954       // ~chi0_i -> chi0_j + Z
 955       fac = kinFac2 * (norm(coupSUSYPtr->OLpp[iNeut1][iNeut2])
 956           + norm(coupSUSYPtr->ORpp[iNeut1][iNeut2]));
 957       fac -= 12.0 * mHat * mf1 * pow2(mf2)
 958            * real(coupSUSYPtr->OLpp[iNeut1][iNeut2]
 959            * conj(coupSUSYPtr->ORpp[iNeut1][iNeut2]));
 960       fac /= pow2(mf2) * (1.0 - s2W);
 961     }
 962     else if(iChar1>0 && id2Abs==24){
 963       // ~chi0_i -> chi+_j + W- (or c.c.)
 964
 965       fac = kinFac2 * (norm(coupSUSYPtr->OL[iNeut1][iChar1])
 966           + norm(coupSUSYPtr->OR[iNeut1][iChar1]));
 967       fac -= 12.0 * mHat * mf1 * pow2(mf2)
 968            * real(coupSUSYPtr->OL[iNeut1][iChar1]
 969            * conj(coupSUSYPtr->OR[iNeut1][iChar1]));
 970       fac /= pow2(mf2);
 971     }
 972     else if(id1Abs > 1000000 && id1Abs%100 < 7 && id2Abs < 7){
 973       // ~chi0_k -> ~q + q
 974       bool idown = (id1Abs%2 == 1);
 975       int iq = (id2Abs + 1 )/ 2;
 976       int isq = (abs(id1Abs)/1000000 == 2) ? (abs(id1Abs)%10+1)/2 + 3
 977                                            : (abs(id1Abs)%10+1)/2;
 978
 979       if(idown){
 980         fac  = kinFac * (norm(coupSUSYPtr->LsddX[isq][iq][iNeut1])
 981              + norm(coupSUSYPtr->RsddX[isq][iq][iNeut1]));
 982         fac += 4.0 * mHat * mf2 * real(coupSUSYPtr->LsddX[isq][iq][iNeut1]
 983              * conj(coupSUSYPtr->RsddX[isq][iq][iNeut1]));
 984       }
 985       else{
 986         fac = kinFac * (norm(coupSUSYPtr->LsuuX[isq][iq][iNeut1])
 987             + norm(coupSUSYPtr->RsuuX[isq][iq][iNeut1]));
 988         fac += 4.0 * mHat * mf2 * real(coupSUSYPtr->LsuuX[isq][iq][iNeut1]
 989              * conj(coupSUSYPtr->RsuuX[isq][iq][iNeut1]));
 990       }
 991       // Extra multiplicative factor of 3 over sleptons
 992       fac *= 6.0/(1 - s2W);
 993     }
 994     else if(id1Abs > 1000000 && id1Abs%100 > 10 && id1Abs%100 < 17
 995       && id2Abs < 17){
 996       // ~chi0_k -> ~l + l
 997       bool idown = id2Abs%2;
 998       int il = (id2Abs - 9)/ 2;
 999       int isl = (abs(id1Abs)/1000000 == 2) ? (abs(id1Abs)%10+1)/2 + 3
1000                                            : (abs(id1Abs)%10+1)/2;
1001
1002       if(idown){
1003         fac  = kinFac * (norm(coupSUSYPtr->LsllX[isl][il][iNeut1])
1004              + norm(coupSUSYPtr->RsllX[isl][il][iNeut1]));
1005         fac += 4.0 * mHat * mf2 * real(coupSUSYPtr->LsllX[isl][il][iNeut1]
1006              * conj(coupSUSYPtr->RsllX[isl][il][iNeut1]));
1007       }
1008       else{
1009         fac = kinFac * (norm(coupSUSYPtr->LsvvX[isl][il][iNeut1]));
1010       }
1011       fac *= 2.0/(1 - s2W);
1012     }
1013     // TODO: Decays in higgs
1014     // Final width for 2-body decays
1015     widNow = fac * preFac * ps ;
1016     if(DEBUG) {
1017       cout<<idRes<<":: id1:"<<id1Abs<<" id2:"<<id2Abs<<" Width: ";
1018       cout<<scientific<<widNow<<endl;
1019     }
1020   }
1021   else {
1022     //RPV 3-body decays
1023     //Case: UDD-type
1024
1025     if(!coupSUSYPtr->isUDD) return;
1026
1027     if(id1Abs < 7 && id2Abs < 7 && id3Abs < 7){
1028
1029       // Check that quarks compatible with UDD structure
1030       if((id1Abs+id2Abs+id3Abs)%2 == 1) return;
1031       double rvfac,m12min,m12max,integral;
1032       int idInt;
1033
1034       // Loop over mass eigenstate in actual propagator
1035       for(int idIntRes = 1; idIntRes <= 6; idIntRes++){
1036         // Right handed field in the UDD coupling
1037         for (int iSq = 1; iSq <= 3; iSq++){
1038           double m1, m2, m3, mixfac1(0.), mixfac2(0.), mixfac3(0.);
1039           int itemp1,itemp2,itemp3;
1040           // up/down-type internal lines
1041           for(int itype = 1; itype<=3; itype++){
1042             //itype = 1: up
1043             //itype = 2: down1
1044             //itype = 3: down2
1045             if(itype ==1 ) idInt = coupSUSYPtr->idSup(idIntRes);
1046             else idInt = coupSUSYPtr->idSdown(idIntRes);
1047             if(id1Abs%2 == 0){
1048               if(itype == 1){
1049                 itemp3 = id1Abs;
1050                 itemp1 = id2Abs;
1051                 itemp2 = id3Abs;
1052                 rvfac = pow2(
1053                   coupSUSYPtr->rvUDD[iSq][(id2Abs+1)/2][(id3Abs+1)/2]);
1054               }else if (itype ==2){
1055                 itemp3 = id2Abs;
1056                 itemp1 = id1Abs;
1057                 itemp2 = id3Abs;
1058                 rvfac = pow2(
1059                   coupSUSYPtr->rvUDD[(id1Abs+1)/2][iSq][(id3Abs+1)/2]);
1060               } else{
1061                 itemp3 = id3Abs;
1062                 itemp1 = id1Abs;
1063                 itemp2 = id2Abs;
1064                 rvfac = pow2(
1065                   coupSUSYPtr->rvUDD[(id1Abs+1)/2][(id2Abs+1)/2][iSq]);
1066               }
1067             }else if(id2Abs%2 == 0){
1068               if(itype==1){
1069                 itemp3 = id2Abs;
1070                 itemp1 = id1Abs;
1071                 itemp2 = id3Abs;
1072                 rvfac = pow2(
1073                   coupSUSYPtr->rvUDD[iSq][(id1Abs+1)/2][(id3Abs+1)/2]);
1074               }else if (itype ==2){
1075                 itemp3 = id1Abs;
1076                 itemp1 = id2Abs;
1077                 itemp2 = id3Abs;
1078                 rvfac = pow2(
1079                   coupSUSYPtr->rvUDD[(id2Abs+1)/2][iSq][(id3Abs+1)/2]);
1080               } else{
1081                 itemp3 = id3Abs;
1082                 itemp1 = id2Abs;
1083                 itemp2 = id1Abs;
1084                 rvfac = pow2(
1085                   coupSUSYPtr->rvUDD[(id2Abs+1)/2][(id2Abs+1)/2][iSq]);
1086               }
1087             }else{
1088               if(itype==1){
1089                 itemp3 = id3Abs;
1090                 itemp1 = id1Abs;
1091                 itemp2 = id2Abs;
1092                 rvfac = pow2(
1093                   coupSUSYPtr->rvUDD[iSq][(id1Abs+1)/2][(id2Abs+1)/2]);
1094               }else if (itype ==2){
1095                 itemp3 = id2Abs;
1096                 itemp1 = id1Abs;
1097                 itemp2 = id3Abs;
1098                 rvfac = pow2(
1099                   coupSUSYPtr->rvUDD[(id3Abs+1)/2][iSq][(id2Abs+1)/2]);
1100               } else{
1101                 itemp3 = id3Abs;
1102                 itemp1 = id1Abs;
1103                 itemp2 = id2Abs;
1104                 rvfac = pow2(
1105                   coupSUSYPtr->rvUDD[(id3Abs+1)/2][(id1Abs+1)/2][iSq]);
1106               }
1107
1108             }
1109
1110             m1 = particleDataPtr->m0(itemp1);
1111             m2 = particleDataPtr->m0(itemp2);
1112             m3 = particleDataPtr->m0(itemp3);
1113
1114             m12min = pow2(m1+m2);
1115             m12max = pow2(mHat-m3);
1116
1117             // Ignore mode when 2-body decay is possible
1118             if(mRes > particleDataPtr->m0(idInt) + particleDataPtr->m0(itemp3))
1119               continue;
1120
1121             // Single diagram squared terms
1122             psi.setInternal(idRes, itemp1, itemp2, itemp3, idInt, 0);
1123             // Mixing with R-states
1124             if(itype == 1)
1125               mixfac1 = norm(coupSUSYPtr->Rusq[idIntRes][iSq+3]);
1126             else
1127               mixfac1 = norm(coupSUSYPtr->Rdsq[idIntRes][iSq+3]);
1128
1129             if(abs(rvfac * mixfac1) > 1.0e-8) {
1130               integral =  integrateGauss(&psi,m12min,m12max,1.0e-4);
1131               widNow += rvfac * mixfac1 * integral;
1132             //   if(DEBUG || idRes == 1000023)
1133             //  cout << scientific << setw(10) <<"Psi: intRes: "<<idInt
1134             //       <<" integral:"<<integral<<" mixfac:"<<mixfac1
1135             //       <<" widNow:"<<widNow<<endl;
1136             }
1137
1138             // Mixing of diagrams with different internal squarks
1139             // of same isospin
1140             for (int idIntRes2 = 1; idIntRes2 <= 6; idIntRes2++){
1141               if(idIntRes2 == idIntRes) continue;
1142               int idInt2;
1143               if(itype == 1 ){
1144                 idInt2 = coupSUSYPtr->idSup(idIntRes2);
1145                 mixfac2 = 2.0 * real(coupSUSYPtr->Rusq[idIntRes][iSq+3]
1146                         * conj(coupSUSYPtr->Rusq[idIntRes2][iSq+3]));
1147               } else {
1148                 idInt2 = coupSUSYPtr->idSdown(idIntRes2);
1149                 mixfac2 = 2.0 * real(coupSUSYPtr->Rdsq[idIntRes][iSq+3]
1150                         * conj(coupSUSYPtr->Rdsq[idIntRes2][iSq+3]));
1151               }
1152
1153               // Ignore mode when 2-body decay is possible
1154               if(mRes > particleDataPtr->m0(idInt2)
1155                       + particleDataPtr->m0(itemp3)) continue;
1156
1157               upsil.setInternal(idRes,itemp1, itemp2,itemp3,idInt,idInt2);
1158               if(abs(rvfac * mixfac2) > 0.0) {
1159                 integral =  integrateGauss(&upsil,m12min,m12max,1.0e-4);
1160                 widNow += rvfac * mixfac2 * integral;
1161                 // if(DEBUG || idRes == 1000023)
1162                 //   cout << scientific << setw(10) <<"Upsilon: intRes: "
1163                 //        <<idInt<<" intRes2:"<<idInt2<<" integral:"<<integral
1164                 //        <<" mixfac:"<<mixfac2<<" widNow:"<<widNow<<endl;
1165               }
1166             }
1167
1168             // Interference between two diagrams with quarks
1169             // of different isospin
1170
1171             for (int idIntRes2 = 1; idIntRes2 <= 6; idIntRes2++){
1172               if(itype != 1 && idIntRes2 == idIntRes) continue;
1173               int idInt2;
1174
1175               for (int iSq2 = 1; iSq2 <= 3; iSq2++){
1176                 if(itype == 1 ){
1177                   idInt2 = coupSUSYPtr->idSdown(idIntRes2);
1178                   mixfac3 = 2.0 * real(coupSUSYPtr->Rusq[idIntRes][iSq+3]
1179                           * conj(coupSUSYPtr->Rdsq[idIntRes2][iSq2+3]));
1180                 } else {
1181                   idInt2 = coupSUSYPtr->idSdown(idIntRes2);
1182                   mixfac3 = 2.0 * real(coupSUSYPtr->Rdsq[idIntRes][iSq+3]
1183                           * conj(coupSUSYPtr->Rdsq[idIntRes2][iSq2+3]));
1184                 }
1185
1186                 if(abs(rvfac * mixfac3) > 0.0) {
1187                   phi.setInternal(idRes,itemp1, itemp2,itemp3,idInt,idInt2);
1188                   //integral = 0;
1189                   //if(idIntRes == 2 && iSq2 ==4)
1190                   integral =  integrateGauss(&phi,m12min,m12max,1.0e-4);
1191                   widNow -= rvfac * mixfac2 * integral;
1192                 }
1193               }
1194             }
1195           }
1196         }
1197       }
1198     }
1199     // Normalisation.  Extra factor of 2 to account for the fact that
1200     // d_i, d_j will always be ordered in ascending order. N_c! = 6
1201     widNow *= 12.0 /(pow3(2.0 * M_PI * mHat) * 32.0);
1202   }
1203
1204   return;
1205 }
1206
1207 //==========================================================================
1208
1209 //  Class ResonanceChar
1210 //  Derived class for Neutralino Resonances
1211 //  Decays into higgses/sleptons not yet implemented
1212
1213 //--------------------------------------------------------------------------
1214
1215 void ResonanceChar::initConstants(){
1216
1217   alpEM = coupSUSYPtr->alphaEM(mHat * mHat);
1218   s2W   = coupSUSYPtr->sin2W;
1219   return;
1220 }
1221
1222 //--------------------------------------------------------------------------
1223
1224 // Calculate various common prefactors for the current mass.
1225 void  ResonanceChar::calcPreFac(bool){
1226
1227   preFac = alpEM / (8.0 * s2W * pow(mHat,3));
1228   return;
1229 }
1230
1231 //--------------------------------------------------------------------------
1232
1233 // Calculate width for currently considered channel.
1234 void  ResonanceChar::calcWidth(bool){
1235
1236   widNow = 0.0;
1237   if(ps == 0.) return;
1238
1239   if(mult ==2){
1240     double fac = 0.0;
1241     kinFac = mHat * mHat - mf1 * mf1 + mf2 * mf2;
1242     kinFac2 = pow(mHat,4) + pow(mf1,4) - 2.0 * pow(mf2,4)
1243       + pow2(mHat) * pow2(mf2) + pow2(mf1)
1244       * pow2(mf2) - 2.0 * pow2(mHat) * pow2(mf1);
1245
1246     int idChar1 = coupSUSYPtr->typeChar(idRes);
1247     int idChar2 = coupSUSYPtr->typeChar(id1Abs);
1248     int idNeut1 = coupSUSYPtr->typeNeut(id1Abs);
1249
1250     if(idChar2>0 && id2Abs == 23){
1251       // ~chi_i -> chi_j + Z
1252       fac = kinFac2 * (norm(coupSUSYPtr->OLp[idChar1][idChar2])
1253           + norm(coupSUSYPtr->ORp[idChar1][idChar2]));
1254       fac -= 12.0 * mHat * mf1 * pow2(mf2)
1255            * real(coupSUSYPtr->OLp[idChar1][idChar2]
1256            * conj(coupSUSYPtr->ORp[idChar1][idChar2]));
1257       fac /= pow2(mf2) * (1.0 - s2W);
1258     }
1259     else if(idNeut1>0 && id2Abs==24){
1260       // ~chi_i -> chi0_j + W- (or c.c.)
1261
1262       fac  = kinFac2 * (norm(coupSUSYPtr->OL[idNeut1][idChar1])
1263            + norm(coupSUSYPtr->OR[idNeut1][idChar1]));
1264       fac -= 12.0 * mHat * mf1 * pow2(mf2)
1265            * real(coupSUSYPtr->OL[idNeut1][idChar1]
1266            * conj(coupSUSYPtr->OR[idNeut1][idChar1]));
1267       fac /= pow2(mf2);
1268     }
1269     else if(id1Abs > 1000000 && id1Abs%100 < 7 && id2Abs < 7){
1270       // ~chi0_k -> ~q + q
1271       bool idown = (id1Abs%2 == 1);
1272       int iq = (id2Abs + 1 )/ 2;
1273       int isq = (abs(id1Abs)/1000000 == 2) ? (abs(id1Abs)%10+1)/2 + 3
1274                                            : (abs(id1Abs)%10+1)/2;
1275
1276       if(idown){
1277         fac  = kinFac * (norm(coupSUSYPtr->LsduX[isq][iq][idChar1])
1278              + norm(coupSUSYPtr->RsduX[isq][iq][idChar1]));
1279         fac += 4.0 * mHat * mf2
1280              * real(coupSUSYPtr->LsduX[isq][iq][idChar1]
1281              * conj(coupSUSYPtr->RsduX[isq][iq][idChar1]));
1282       }
1283       else{
1284         fac  = kinFac * (norm(coupSUSYPtr->LsudX[isq][iq][idChar1])
1285              + norm(coupSUSYPtr->RsudX[isq][iq][idChar1]));
1286         fac += 4.0 * mHat * mf2
1287              * real(coupSUSYPtr->LsudX[isq][iq][idChar1]
1288              * conj(coupSUSYPtr->RsudX[isq][iq][idChar1]));
1289       }
1290       fac *= 6.0/(1 - s2W);
1291     }
1292     else if(id1Abs > 1000000 && id1Abs%100 > 10 && id1Abs%100 < 17
1293       && id2Abs < 17){
1294       // ~chi+_k -> ~l + l
1295       bool idown = id2Abs%2;
1296       int il = (id2Abs - 9)/ 2;
1297       int isl = (abs(id1Abs)/1000000 == 2) ? (abs(id1Abs)%10+1)/2 + 3
1298                                            : (abs(id1Abs)%10+1)/2;
1299
1300       if(idown){
1301         fac  = kinFac * (norm(coupSUSYPtr->LslvX[isl][il][idChar1])
1302              + norm(coupSUSYPtr->RslvX[isl][il][idChar1]));
1303         fac += 4.0 * mHat * mf2 * real(coupSUSYPtr->LslvX[isl][il][idChar1]
1304              * conj(coupSUSYPtr->RslvX[isl][il][idChar1]));
1305       }
1306       else{
1307         fac = kinFac * (norm(coupSUSYPtr->LsvlX[isl][il][idChar1]));
1308       }
1309       fac *= 2.0/(1 - s2W);
1310     }
1311
1312     // TODO: Decays in higgs
1313     widNow = fac * preFac * ps ;
1314     if(DEBUG) {
1315       cout<<idRes<<":: id1:"<<id1Abs<<" id2:"<<id2Abs<<" Width: ";
1316       cout<<scientific<<widNow<<endl;
1317     }
1318   }else{
1319     //TODO: Implement Chargino 3-body decays
1320   }
1321   return;
1322 }
1323
1324
1325 //==========================================================================
1326 // The ResonanceSlepton class
1327 // Derived class for Slepton (and sneutrino) resonances
1328
1329 //--------------------------------------------------------------------------
1330
1331 // Initialize constants.
1332
1333 void ResonanceSlepton::initConstants() {
1334
1335   // Locally stored properties and couplings.
1336   alpEM = coupSUSYPtr->alphaEM(mHat * mHat);
1337   s2W   = coupSUSYPtr->sin2W;
1338 }
1339
1340 //--------------------------------------------------------------------------
1341
1342 // Calculate various common prefactors for the current mass.
1343
1344 void ResonanceSlepton::calcPreFac(bool) {
1345
1346   // Common coupling factors.
1347   preFac = 1.0 / (s2W * pow(mHat,3));
1348
1349 }
1350
1351 //--------------------------------------------------------------------------
1352
1353 // Calculate width for currently considered channel.
1354
1355 void ResonanceSlepton::calcWidth(bool) {
1356
1357   // Slepton type -- in u_i/d_i and generation
1358   int ksusy = 1000000;
1359   int isl = (abs(idRes)/ksusy == 2) ? (abs(idRes)%10+1)/2 + 3
1360                                     : (abs(idRes)%10+1)/2;
1361   int il = (id2Abs-9)/2;
1362   bool islep = idRes%2;
1363
1364   // Check that mass is above threshold.
1365   if(ps == 0.) return;
1366   else{
1367     // Two-body decays
1368     kinFac = (mHat * mHat - mf1 * mf1 - mf2 * mf2);
1369
1370     double fac = 0.0 , wid = 0.0;
1371
1372     //Case 1: RPV: To be implemented
1373     //Case 2: slepton + gaugino
1374
1375     if (id1Abs > ksusy && id2Abs > 10 && id2Abs < 17) {
1376       for(int i=1; i<6 ; i++){
1377         // ~ell/~nu -> ~chi0 + ell/nu
1378         if(coupSUSYPtr->idNeut(i)==id1Abs && idRes%2 == id2Abs%2){
1379           fac = alpEM *  preFac / (2.0 * (1 - s2W));
1380           if(islep)
1381             wid = kinFac * (norm(coupSUSYPtr->LsllX[isl][il][i])
1382                 + norm(coupSUSYPtr->RsllX[isl][il][i]))
1383                 - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsllX[isl][il][i]
1384                 * conj(coupSUSYPtr->RsllX[isl][il][i]));
1385           else
1386             wid = kinFac * (norm(coupSUSYPtr->LsvvX[isl][il][i])
1387                 + norm(coupSUSYPtr->RsvvX[isl][il][i]))
1388                 - 4.0 * mHat * mf2 * real(coupSUSYPtr->LsvvX[isl][il][i]
1389                 * conj(coupSUSYPtr->RsvvX[isl][il][i]));
1390         }
1391
1392         // ~ell/~nu -> ~chi- + nu/ell
1393         else if (i < 3 && coupSUSYPtr->idChar(i)==id1Abs
1394           && idRes%2 != id2Abs%2){
1395
1396           fac = alpEM *  preFac / (4.0 * (1 - s2W));
1397           if(islep)
1398             wid = kinFac * (norm(coupSUSYPtr->LslvX[isl][il][i])
1399                 + norm(coupSUSYPtr->RslvX[isl][il][i]))
1400                 - 4.0 * mHat * mf2 * real(coupSUSYPtr->LslvX[isl][il][i]
1401                 * conj(coupSUSYPtr->RslvX[isl][il][i]));
1402           else
1403             wid = kinFac * (norm(coupSUSYPtr->LslvX[isl][il][i])
1404                 + norm(coupSUSYPtr->RslvX[isl][il][i]))
1405                 - 4.0 * mHat * mf2 * real(coupSUSYPtr->LslvX[isl][il][i]
1406                 * conj(coupSUSYPtr->RslvX[isl][il][i]));
1407         }
1408       }
1409     }
1410
1411     //Case 3: ~l_i -> ~l_j + Z/W
1412     else if (id1Abs > ksusy+10 && id1Abs%100 < 17
1413       && (id2Abs == 23 || id2Abs == 24)){
1414
1415       // factor of lambda^(3/2) = ps^3;
1416       fac = alpEM * preFac/(16.0 * pow2(mf2) * (1.0 - s2W)) * pow2(ps) ;
1417
1418       int isl2 = (id1Abs/ksusy == 2) ? (id1Abs%10+1)/2 + 3: (id1Abs%10+1)/2;
1419
1420       if(id2Abs == 23 && id1Abs%2 == idRes%2){
1421         if(islep)
1422           wid = norm(coupSUSYPtr->LslslZ[isl][isl2]
1423               + coupSUSYPtr->RslslZ[isl][isl2]);
1424         else
1425           wid = norm(coupSUSYPtr->LsvsvZ[isl][isl2]
1426               + coupSUSYPtr->RsvsvZ[isl][isl2]);
1427       }
1428       else if (id2Abs == 24 && id1Abs%2 != idRes%2){
1429         if(islep)
1430           wid = norm(coupSUSYPtr->LslsvW[isl2][isl]);
1431         else
1432           wid = norm(coupSUSYPtr->LslsvW[isl][isl2]);
1433       }
1434     }
1435
1436     // TODO: Case ~l_i -> ~l_j + h/H
1437
1438
1439     widNow = fac * wid * ps;
1440     if(DEBUG) cout<<idRes<<":: id1:"<<id1Abs<<" id2:"<<id2Abs
1441                   <<" Width: "<<widNow<<endl;
1442     return;
1443   }
1444
1445 }
1446
1447 //==========================================================================
1448
1449 // Gaussian Integrator for 3-body decay widths
1450
1451 double SUSYResonanceWidths::integrateGauss(WidthFunction* widthFn,
1452   double xlo, double xhi, double tol) {
1453
1454   //8-point unweighted
1455   static double x8[4]={0.96028985649753623,
1456                        0.79666647741362674,
1457                        0.52553240991632899,
1458                        0.18343464249564980};
1459   static double w8[4]={0.10122853629037626,
1460                        0.22238103445337447,
1461                        0.31370664587788729,
1462                        0.36268378337836198};
1463   //16-point unweighted
1464   static double x16[8]={0.98940093499164993,
1465                        0.94457502307323258,
1466                        0.86563120238783174,
1467                        0.75540440835500303,
1468                        0.61787624440264375,
1469                        0.45801677765722739,
1470                        0.28160355077925891,
1471                        0.09501250983763744};
1472   static double w16[8]={0.027152459411754095,
1473                        0.062253523938647893,
1474                        0.095158511682492785,
1475                        0.12462897125553387,
1476                        0.14959598881657673,
1477                        0.16915651939500254,
1478                        0.18260341504492359,
1479                        0.18945061045506850};
1480   //boundary checks
1481   if (xlo == xhi) {
1482     cerr<<"xlo = xhi"<<endl;
1483     return 0.0;
1484   }
1485   if ( xlo > xhi ) {
1486     cerr<<" (integrateGauss:) -> xhi < xlo"<<endl;
1487     return 0.0;
1488   }
1489   //initialize
1490   double sum = 0.0;
1491   double c = 0.001/abs(xhi-xlo);
1492   double zlo = xlo;
1493   double zhi = xhi;
1494
1495   bool nextbin = true;
1496
1497   while ( nextbin ) {
1498
1499     double zmi = 0.5*(zhi+zlo); // midpoint
1500     double zmr = 0.5*(zhi-zlo); // midpoint, relative to zlo
1501
1502     //calculate 8-point and 16-point quadratures
1503     double s8=0.0;
1504     for (int i=0;i<4;i++) {
1505       double dz = zmr * x8[i];
1506       s8 += w8[i]*(widthFn->function(zmi+dz) + widthFn->function(zmi-dz));
1507     }
1508     s8 *= zmr;
1509     double s16=0.0;
1510     for (int i=0;i<8;i++) {
1511       double dz = zmr * x16[i];
1512       s16 += w16[i]*(widthFn->function(zmi+dz) + widthFn->function(zmi-dz));
1513     }
1514     s16 *= zmr;
1515     if (abs(s16-s8) < tol*(1+abs(s16))) {
1516       //precision in this bin OK, add to cumulative and go to next
1517       nextbin=true;
1518       sum += s16;
1519       //next bin: LO = end of current, HI = end of integration region.
1520       zlo=zhi;
1521       zhi=xhi;
1522       if ( zlo == zhi ) nextbin = false;
1523     } else {
1524       //precision in this bin not OK, subdivide.
1525       if (1.0 + c*abs(zmr) == 1.0) {
1526         cerr << " (integrateGauss:) too high accuracy required"<<endl;
1527         sum = 0.0 ;
1528         break;
1529       }
1530       zhi=zmi;
1531       nextbin=true;
1532     }
1533   }
1534
1535   return sum;
1536 }
1537
1538 //======================================================================
1539
1540 } //end namespace Pythia8
1541
1542