PYTHIA8/pythia8145/src/SigmaSUSY.cxx

   1 // SigmaSUSY.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 // supersymmetry simulation classes.
   8
   9 #include "SigmaSUSY.h"
  10
  11 namespace Pythia8 {
  12
  13 //==========================================================================
  14
  15 // Sigma2qqbar2chi0chi0
  16 // Cross section for gaugino pair production: neutralino pair
  17
  18 //--------------------------------------------------------------------------
  19
  20 // Initialize process.
  21
  22 void Sigma2qqbar2chi0chi0::initProc() {
  23
  24   //Typecast to the correct couplings
  25   coupSUSYPtr = (CoupSUSY*) couplingsPtr;
  26
  27   // Construct name of process.
  28   nameSave = "q qbar' -> " + particleDataPtr->name(id3) + " "
  29     + particleDataPtr->name(id4);
  30
  31   // Secondary open width fraction.
  32   openFracPair = particleDataPtr->resOpenFrac(id3, id4);
  33
  34 }
  35
  36 //--------------------------------------------------------------------------
  37
  38 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
  39
  40 void Sigma2qqbar2chi0chi0::sigmaKin() {
  41
  42   // Common flavour-independent factor.
  43   sigma0 = M_PI /3.0/ sH2 / pow2(coupSUSYPtr->sin2W) * pow2(alpEM)
  44     * openFracPair;
  45
  46   // Auxiliary factors for use below
  47   ui       = uH - s3;
  48   uj       = uH - s4;
  49   ti       = tH - s3;
  50   tj       = tH - s4;
  51   double sV= sH - pow2(coupSUSYPtr->mZpole);
  52   double d = pow2(sV) + pow2(coupSUSYPtr->mZpole * coupSUSYPtr->wZpole);
  53   propZ    = complex( sV / d, coupSUSYPtr->mZpole * coupSUSYPtr->wZpole / d);
  54
  55 }
  56
  57 //--------------------------------------------------------------------------
  58
  59 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
  60
  61 double Sigma2qqbar2chi0chi0::sigmaHat() {
  62
  63   // Only allow quark-antiquark incoming states
  64   if (id1*id2 >= 0) {
  65     return 0.0;
  66   }
  67
  68   // Only allow incoming states with sum(charge) = 0
  69   if ((id1+id2) % 2 != 0) {
  70     return 0.0;
  71   }
  72
  73   if(id1<0) swapTU = true;
  74
  75   // Shorthands
  76   int idAbs1    = abs(id1);
  77   int idAbs2    = abs(id2);
  78
  79   // Flavour-dependent kinematics-dependent couplings.
  80   complex QuLL(0.0),QtLL(0.0),QuRR(0.0),QtRR(0.0);
  81   complex QuLR(0.0),QtLR(0.0),QuRL(0.0),QtRL(0.0);
  82
  83   // s-channel Z couplings
  84   if (idAbs1 == idAbs2) {
  85     QuLL = coupSUSYPtr->LqqZ[idAbs1] * coupSUSYPtr->OLpp[id3chi][id4chi] * propZ / 2.0;
  86     QtLL = coupSUSYPtr->LqqZ[idAbs1] * coupSUSYPtr->ORpp[id3chi][id4chi] * propZ / 2.0;
  87     QuRR = coupSUSYPtr->RqqZ[idAbs1] * coupSUSYPtr->ORpp[id3chi][id4chi] * propZ / 2.0;
  88     QtRR = coupSUSYPtr->RqqZ[idAbs1] * coupSUSYPtr->OLpp[id3chi][id4chi] * propZ / 2.0;
  89   }
  90
  91   // Flavour indices
  92   int ifl1 = (idAbs1+1) / 2;
  93   int ifl2 = (idAbs2+1) / 2;
  94
  95   // Add t-channel squark flavour sums to QmXY couplings
  96   for (int ksq=1; ksq<=6; ksq++) {
  97
  98     // squark id and squark-subtracted u and t
  99     int idsq=((ksq+2)/3)*1000000 + 2*((ksq-1) % 3) + (idAbs1+1) % 2 + 1;
 100     double msq2    = pow(particleDataPtr->m0(idsq),2);
 101     double usq     = uH - msq2;
 102     double tsq     = tH - msq2;
 103
 104     // Couplings
 105     complex Lsqq1X3 = coupSUSYPtr->LsuuX[ksq][ifl1][id3chi];
 106     complex Lsqq1X4 = coupSUSYPtr->LsuuX[ksq][ifl1][id4chi];
 107     complex Lsqq2X3 = coupSUSYPtr->LsuuX[ksq][ifl2][id3chi];
 108     complex Lsqq2X4 = coupSUSYPtr->LsuuX[ksq][ifl2][id4chi];
 109     complex Rsqq1X3 = coupSUSYPtr->RsuuX[ksq][ifl1][id3chi];
 110     complex Rsqq1X4 = coupSUSYPtr->RsuuX[ksq][ifl1][id4chi];
 111     complex Rsqq2X3 = coupSUSYPtr->RsuuX[ksq][ifl2][id3chi];
 112     complex Rsqq2X4 = coupSUSYPtr->RsuuX[ksq][ifl2][id4chi];
 113     if (idAbs1 % 2 != 0) {
 114       Lsqq1X3 = coupSUSYPtr->LsddX[ksq][ifl1][id3chi];
 115       Lsqq1X4 = coupSUSYPtr->LsddX[ksq][ifl1][id4chi];
 116       Lsqq2X3 = coupSUSYPtr->LsddX[ksq][ifl2][id3chi];
 117       Lsqq2X4 = coupSUSYPtr->LsddX[ksq][ifl2][id4chi];
 118       Rsqq1X3 = coupSUSYPtr->RsddX[ksq][ifl1][id3chi];
 119       Rsqq1X4 = coupSUSYPtr->RsddX[ksq][ifl1][id4chi];
 120       Rsqq2X3 = coupSUSYPtr->RsddX[ksq][ifl2][id3chi];
 121       Rsqq2X4 = coupSUSYPtr->RsddX[ksq][ifl2][id4chi];
 122     }
 123
 124     // QuXY
 125     QuLL += conj(Lsqq1X4)*Lsqq2X3/usq;
 126     QuRR += conj(Rsqq1X4)*Rsqq2X3/usq;
 127     QuLR += conj(Lsqq1X4)*Rsqq2X3/usq;
 128     QuRL += conj(Rsqq1X4)*Lsqq2X3/usq;
 129
 130
 131     // QtXY
 132     QtLL -= conj(Lsqq1X3)*Lsqq2X4/tsq;
 133     QtRR -= conj(Rsqq1X3)*Rsqq2X4/tsq;
 134     QtLR += conj(Lsqq1X3)*Rsqq2X4/tsq;
 135     QtRL += conj(Rsqq1X3)*Lsqq2X4/tsq;
 136
 137   }
 138
 139   // Overall factor multiplying each coupling; multiplied at the end as fac^2
 140   double fac = (1.0-coupSUSYPtr->sin2W);
 141   if(abs(id3)==abs(id4)) fac *= sqrt(2); // for identical final state particles
 142
 143   // Compute matrix element weight
 144   double weight = 0;
 145   double facLR = uH*tH - s3*s4;
 146   double facMS = m3*m4*sH;
 147
 148   // Average over separate helicity contributions
 149   // LL (ha = -1, hb = +1) (divided by 4 for average)
 150   weight += norm(QuLL) * ui * uj + norm(QtLL) * ti * tj
 151     + 2 * real(conj(QuLL) * QtLL) * facMS;
 152   // RR (ha =  1, hb = -1) (divided by 4 for average)
 153   weight += norm(QtRR) * ti * tj + norm(QuRR) * ui * uj
 154     + 2 * real(conj(QuRR) * QtRR) * facMS;
 155   // RL (ha =  1, hb =  1) (divided by 4 for average)
 156   weight += norm(QuRL) * ui * uj + norm(QtRL) * ti * tj
 157     + real(conj(QuRL) * QtRL) * facLR;
 158   // LR (ha = -1, hb = -1) (divided by 4 for average)
 159   weight += norm(QuLR) * ui * uj + norm(QtLR) * ti * tj
 160     + real(conj(QuLR) * QtLR) * facLR;
 161
 162   // Cross section, including colour factor.
 163   double sigma = sigma0 * weight / pow2(fac);
 164
 165   // Answer.
 166   return sigma;
 167
 168 }
 169
 170 //--------------------------------------------------------------------------
 171
 172 // Select identity, colour and anticolour.
 173
 174 void Sigma2qqbar2chi0chi0::setIdColAcol() {
 175
 176   // Set flavours.
 177   setId( id1, id2, id3, id4);
 178
 179   // Colour flow topologies. Swap when antiquarks.
 180   if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
 181   else              setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
 182   if (id1 < 0) swapColAcol();
 183
 184 }
 185
 186 //==========================================================================
 187
 188 // Sigma2qqbar2charchi0
 189 // Cross section for gaugino pair production: neutralino-chargino
 190
 191 //--------------------------------------------------------------------------
 192
 193 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
 194
 195 void Sigma2qqbar2charchi0::sigmaKin() {
 196
 197   // Common flavour-independent factor.
 198
 199   sigma0 = M_PI / sH2 / 3.0 / pow2(coupSUSYPtr->sin2W) * pow2(alpEM) ;
 200   sigma0 /= 2.0 * (1 - coupSUSYPtr->sin2W) ;
 201
 202   // Auxiliary factors for use below
 203   ui        = uH - s3;
 204   uj        = uH - s4;
 205   ti        = tH - s3;
 206   tj        = tH - s4;
 207   double sW = sH - pow2(coupSUSYPtr->mWpole);
 208   double d  = pow2(sW) + pow2(coupSUSYPtr->mWpole * coupSUSYPtr->wWpole);
 209   propW     = complex( sW / d, coupSUSYPtr->mWpole * coupSUSYPtr->wWpole / d);
 210
 211 }
 212
 213 //--------------------------------------------------------------------------
 214
 215 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
 216
 217 double Sigma2qqbar2charchi0::sigmaHat() {
 218
 219   // Only allow particle-antiparticle incoming states
 220   if (id1*id2 >= 0) {
 221     return 0.0;
 222   }
 223
 224   // Only allow incoming states with sum(charge) = final state
 225   if (abs(id1) % 2 == abs(id2) % 2) return 0.0;
 226   int isPos  = (id3chi > 0 ? 1 : 0);
 227   if (id1 < 0 && id1 > -10 && abs(id1) % 2 == 1-isPos ) return 0.0;
 228   else if (id1 > 0 && id1 < 10 && abs(id1) % 2 == isPos ) return 0.0;
 229
 230   // Flavour-dependent kinematics-dependent couplings.
 231   int idAbs1  = abs(id1);
 232   int iChar = abs(id3chi);
 233   int iNeut = abs(id4chi);
 234
 235   complex QuLL(0.0),QtLL(0.0),QuRR(0.0),QtRR(0.0);
 236   complex QuLR(0.0),QtLR(0.0),QuRL(0.0),QtRL(0.0);
 237
 238   // Calculate everything from udbar -> ~chi+ ~chi0 template process
 239
 240   // u dbar , ubar d : do nothing
 241   // dbar u , d ubar : swap 1<->2 and t<->u
 242
 243   int iGu = abs(id1)/2;
 244   int iGd = abs(id2+1)/2;
 245
 246   if (idAbs1 % 2 != 0) {
 247     swapTU = true;
 248     iGu = abs(id2)/2;
 249     iGd = abs(id1+1)/2;
 250   }
 251
 252   // s-channel W contribution
 253   QuLL = conj(coupSUSYPtr->LudW[iGu][iGd])
 254     * conj(coupSUSYPtr->OL[iNeut][abs(iChar)])
 255     * propW / sqrt(2.0);
 256   QtLL = conj(coupSUSYPtr->LudW[iGu][iGd])
 257     * conj(coupSUSYPtr->OR[iNeut][iChar])
 258     * propW / sqrt(2.0);
 259
 260   // Add t-channel squark flavour sums to QmXY couplings
 261   for (int jsq=1; jsq<=6; jsq++) {
 262
 263     int idsu=((jsq+2)/3)*1000000 + 2*((jsq-1) % 3) + 2;
 264     int idsd=((jsq+2)/3)*1000000 + 2*((jsq-1) % 3) + 1;
 265     double msd2 = pow(particleDataPtr->m0(idsd),2);
 266     double msu2 = pow(particleDataPtr->m0(idsu),2);
 267     double tsq  = tH - msd2;
 268     double usq  = uH - msu2;
 269
 270     QuLL += conj(coupSUSYPtr->LsuuX[jsq][iGu][iNeut])
 271       * conj(coupSUSYPtr->LsudX[jsq][iGd][iChar])/usq;
 272     QuLR += conj(coupSUSYPtr->LsuuX[jsq][iGu][iNeut])
 273       * conj(coupSUSYPtr->RsudX[jsq][iGd][iChar])/usq;
 274     QuRR += conj(coupSUSYPtr->RsuuX[jsq][iGu][iNeut])
 275       * conj(coupSUSYPtr->RsudX[jsq][iGd][iChar])/usq;
 276     QuRL += conj(coupSUSYPtr->RsuuX[jsq][iGu][iNeut])
 277       * conj(coupSUSYPtr->LsudX[jsq][iGd][iChar])/usq;
 278
 279     QtLL -= conj(coupSUSYPtr->LsduX[jsq][iGu][iChar])
 280       * coupSUSYPtr->LsddX[jsq][iGd][iNeut]/tsq;
 281     QtRR -= conj(coupSUSYPtr->RsduX[jsq][iGu][iChar])
 282       * coupSUSYPtr->RsddX[jsq][iGd][iNeut]/tsq;
 283     QtLR += conj(coupSUSYPtr->LsduX[jsq][iGu][iChar])
 284       * coupSUSYPtr->RsddX[jsq][iGd][iNeut]/tsq;
 285     QtRL += conj(coupSUSYPtr->RsduX[jsq][iGu][iChar])
 286       * coupSUSYPtr->LsddX[jsq][iGd][iNeut]/tsq;
 287   }
 288
 289   // Compute matrix element weight
 290   double weight = 0;
 291
 292   // Average over separate helicity contributions
 293   // (if swapped, swap ha, hb if computing polarized cross sections)
 294   // LL (ha = -1, hb = +1) (divided by 4 for average)
 295   weight += norm(QuLL) * ui * uj + norm(QtLL) * ti * tj
 296     + 2 * real(conj(QuLL) * QtLL) * m3 * m4 * sH;
 297   // RR (ha =  1, hb = -1) (divided by 4 for average)
 298   weight += norm(QtRR) * ti * tj + norm(QuRR) * ui * uj
 299     + 2 * real(conj(QuRR) * QtRR) * m3 * m4 * sH;
 300   // RL (ha =  1, hb =  1) (divided by 4 for average)
 301   weight += norm(QuRL) * ui * uj + norm(QtRL) * ti * tj
 302     + real(conj(QuRL) * QtRL) * (uH * tH - s3 * s4);
 303   // LR (ha = -1, hb = -1) (divided by 4 for average)
 304   weight += norm(QuLR) * ui * uj + norm(QtLR) * ti * tj
 305     + real(conj(QuLR) * QtLR) * (uH * tH - s3 * s4);
 306
 307   // Cross section, including colour factor.
 308   double sigma = sigma0 * weight;
 309   if (idAbs1 < 9) sigma /= 3.;
 310
 311   // Answer.
 312   return sigma;
 313
 314 }
 315
 316 //==========================================================================
 317
 318 // Sigma2qqbar2charchar
 319 // Cross section for gaugino pair production: chargino-chargino
 320
 321 //--------------------------------------------------------------------------
 322
 323 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
 324
 325 void Sigma2qqbar2charchar::sigmaKin() {
 326
 327   // Common flavour-independent factor.
 328   sigma0 = M_PI / 3.0 / sH2 / pow2(coupSUSYPtr->sin2W) * pow2(alpEM) ;
 329
 330   // Auxiliary factors for use below
 331   ui       = uH - s3;
 332   uj       = uH - s4;
 333   ti       = tH - s3;
 334   tj       = tH - s4;
 335   double sV= sH - pow2(coupSUSYPtr->mZpole);
 336   double d = pow2(sV) + pow2(coupSUSYPtr->mZpole * coupSUSYPtr->wZpole);
 337   propZ    = complex( sV / d, coupSUSYPtr->mZpole * coupSUSYPtr->wZpole / d);
 338
 339 }
 340
 341 //--------------------------------------------------------------------------
 342
 343 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
 344
 345 double Sigma2qqbar2charchar::sigmaHat() {
 346
 347   // Only allow quark-antiquark incoming states
 348   if (id1*id2 >= 0) {
 349     return 0.0;
 350   }
 351
 352   // Only allow incoming states with sum(charge) = 0
 353   if ((id1+id2) % 2 != 0) {
 354     return 0.0;
 355   }
 356
 357   //if (id1 > 0 || id1==-1 || id1==-3 || id1==-5) return 0.0;
 358   //if (id1 < 0 || id1==1 || id1==3 || id1==5) return 0.0;
 359
 360   swapTU = (id1 < 0 ? true : false);
 361
 362   // Flavour-dependent kinematics-dependent couplings.
 363   int idAbs1    = abs(id1);
 364   int idAbs2    = abs(id2);
 365   int i3        = abs(id3chi);
 366   int i4        = abs(id4chi);
 367
 368   // Flavour-dependent kinematics-dependent couplings.
 369   complex QuLL(0.0),QtLL(0.0),QuRR(0.0),QtRR(0.0);
 370   complex QuLR(0.0),QtLR(0.0),QuRL(0.0),QtRL(0.0);
 371
 372   // Add Z/gamma* for same-flavour in-quarks
 373   if (idAbs1 == idAbs2) {
 374
 375     QuLL = -coupSUSYPtr->LqqZ[idAbs1]*conj(coupSUSYPtr->ORp[i3][i4]);
 376     QtLL = -coupSUSYPtr->LqqZ[idAbs1]*conj(coupSUSYPtr->OLp[i3][i4]);
 377     QuRR = -coupSUSYPtr->RqqZ[idAbs1]*conj(coupSUSYPtr->OLp[i3][i4]);
 378     QtRR = -coupSUSYPtr->RqqZ[idAbs1]*conj(coupSUSYPtr->ORp[i3][i4]);
 379
 380     QuLL *= propZ / 2.0 / (1.0-coupSUSYPtr->sin2W);
 381     QtLL *= propZ / 2.0 / (1.0-coupSUSYPtr->sin2W);
 382     QuRR *= propZ / 2.0 / (1.0-coupSUSYPtr->sin2W);
 383     QtRR *= propZ / 2.0 / (1.0-coupSUSYPtr->sin2W);
 384
 385     // s-channel gamma* (only for same-type charginos)
 386     if (i3 == i4) {
 387
 388       // Charge of in-particles
 389       double q = 2.0/3.0;
 390       if (idAbs1 % 2 == 1) q = -1.0/3.0;
 391       QuLL += q * coupSUSYPtr->sin2W / sH;
 392       QuRR += q * coupSUSYPtr->sin2W / sH;
 393       QtLL += q * coupSUSYPtr->sin2W / sH;
 394       QtRR += q * coupSUSYPtr->sin2W / sH;
 395
 396     }
 397   }
 398
 399   int iG1    = (abs(id1)+1)/2;
 400   int iG2    = (abs(id2)+1)/2;
 401
 402   // Add t- or u-channel squark flavour sums to QmXY couplings
 403   for (int ksq=1; ksq<=6; ksq++) {
 404
 405     if(id1 % 2 == 0) {
 406
 407       // u-channel diagrams only
 408       // up-type incoming; u-channel ~d
 409
 410       int idsd    = ((ksq+2)/3)*1000000 + 2*((ksq-1) % 3) + 1;
 411       double msq  = particleDataPtr->m0(idsd);
 412       double ufac = 2.0 * (uH - pow2(msq));
 413
 414       //u-ubar -> chi-chi+
 415       QuLL += coupSUSYPtr->LsduX[ksq][iG2][i3] * conj(coupSUSYPtr->LsduX[ksq][iG1][i4]) / ufac;
 416       QuRR += coupSUSYPtr->RsduX[ksq][iG2][i3] * conj(coupSUSYPtr->RsduX[ksq][iG1][i4]) / ufac;
 417       QuLR += coupSUSYPtr->RsduX[ksq][iG2][i3] * conj(coupSUSYPtr->LsduX[ksq][iG1][i4]) / ufac;
 418       QuRL += coupSUSYPtr->LsduX[ksq][iG2][i3] * conj(coupSUSYPtr->RsduX[ksq][iG1][i4]) / ufac;
 419
 420
 421     }else{
 422       // t-channel diagrams only;
 423       // down-type incoming; t-channel ~u
 424
 425       int idsu    = ((ksq+2)/3)*1000000 + 2*((ksq-1) % 3) + 2;
 426       double msq  = particleDataPtr->m0(idsu);
 427       double tfac = 2.0 * (tH - pow2(msq));
 428
 429       //d-dbar -> chi-chi+
 430       QtLL -= coupSUSYPtr->LsudX[ksq][iG1][i3] * conj(coupSUSYPtr->LsudX[ksq][iG2][i4]) / tfac;
 431       QtRR -= coupSUSYPtr->RsudX[ksq][iG1][i3] * conj(coupSUSYPtr->RsudX[ksq][iG2][i4]) / tfac;
 432       QtLR += coupSUSYPtr->LsudX[ksq][iG1][i3] * conj(coupSUSYPtr->RsudX[ksq][iG2][i4]) / tfac;
 433       QtRL += coupSUSYPtr->RsudX[ksq][iG1][i3] * conj(coupSUSYPtr->LsudX[ksq][iG2][i4]) / tfac;
 434
 435     }
 436   }
 437    // Compute matrix element weight
 438    double weight = 0;
 439
 440   // Average over separate helicity contributions
 441   // LL (ha = -1, hb = +1) (divided by 4 for average)
 442   weight += norm(QuLL) * ui * uj + norm(QtLL) * ti * tj
 443     + 2 * real(conj(QuLL) * QtLL) * m3 * m4 * sH;
 444   // RR (ha =  1, hb = -1) (divided by 4 for average)
 445   weight += norm(QtRR) * ti * tj + norm(QuRR) * ui * uj
 446     + 2 * real(conj(QuRR) * QtRR) * m3 * m4 * sH;
 447   // RL (ha =  1, hb =  1) (divided by 4 for average)
 448   weight += norm(QuRL) * ui * uj + norm(QtRL) * ti * tj
 449     + real(conj(QuRL) * QtRL) * (uH * tH - s3 * s4);
 450   // LR (ha = -1, hb = -1) (divided by 4 for average)
 451   weight += norm(QuLR) * ui * uj + norm(QtLR) * ti * tj
 452     + real(conj(QuLR) * QtLR) * (uH * tH - s3 * s4);
 453
 454   // Cross section, including colour factor.
 455   double sigma = sigma0 * weight;
 456
 457   // Answer.
 458   return sigma;
 459
 460 }
 461
 462 //==========================================================================
 463
 464 // Sigma2qgchi0squark
 465 // Cross section for gaugino-squark production: neutralino-squark
 466
 467 //--------------------------------------------------------------------------
 468
 469 // Initialize process.
 470
 471 void Sigma2qg2chi0squark::initProc() {
 472
 473   //Typecast to the correct couplings
 474   coupSUSYPtr = (CoupSUSY*) couplingsPtr;
 475
 476   // Construct name of process.
 477   if (id4 % 2 == 0) {
 478     nameSave = "q g -> " + particleDataPtr->name(id3) + " "
 479       + particleDataPtr->name(id4) + " + c.c. (q=u,c)";
 480   }
 481   else {
 482     nameSave = "q g -> " + particleDataPtr->name(id3) + " "
 483       + particleDataPtr->name(id4) + " + c.c. (q=d,s,b)";
 484   }
 485
 486   // Secondary open width fraction.
 487   openFracPair = particleDataPtr->resOpenFrac(id3, id4);
 488
 489 }
 490
 491 //--------------------------------------------------------------------------
 492
 493 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
 494
 495 void Sigma2qg2chi0squark::sigmaKin() {
 496
 497   // Common flavour-independent factor.
 498   // tmp: alphaS = 0.1 for counter-checks
 499   sigma0 = M_PI / sH2 / coupSUSYPtr->sin2W * alpEM * alpS
 500     * openFracPair;
 501
 502   // Auxiliary factors for use below
 503   ui       = uH - s3;
 504   uj       = uH - s4;
 505   ti       = tH - s3;
 506   tj       = tH - s4;
 507
 508 }
 509
 510 //--------------------------------------------------------------------------
 511
 512 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
 513
 514 double Sigma2qg2chi0squark::sigmaHat() {
 515
 516   // Antiquark -> antisquark
 517   int idq = id1;
 518   if (id1 == 21 || id1 == 22) idq = id2;
 519   if (idq < 0) {
 520     id4 = -abs(id4);
 521   } else {
 522     id4 = abs(id4);
 523   }
 524
 525   // tmp: only allow incoming quarks on side 1
 526   //  if (id1 < 0 || id1 == 21) return 0.0;
 527
 528   // Generation index
 529   int iGq = (abs(idq)+1)/2;
 530
 531   // Only accept u(bar) -> ~u(bar) and d(bar) -> ~d(bar)
 532   if (particleDataPtr->chargeType(idq) != particleDataPtr->chargeType(id4))
 533     return 0.0;
 534
 535   // Couplings
 536   complex LsqqX, RsqqX;
 537   if (idq % 2 == 0) {
 538     LsqqX = coupSUSYPtr->LsuuX[id4sq][iGq][id3chi];
 539     RsqqX = coupSUSYPtr->RsuuX[id4sq][iGq][id3chi];
 540   }
 541   else {
 542     LsqqX = coupSUSYPtr->LsddX[id4sq][iGq][id3chi];
 543     RsqqX = coupSUSYPtr->RsddX[id4sq][iGq][id3chi];
 544   }
 545
 546   // Prefactors : swap u and t if gq instead of qg
 547   double fac1, fac2;
 548   if (idq == id1) {
 549     fac1 = -ui/sH + 2.0 * ( uH*tH - s4*s3 )/sH/tj;
 550     fac2 = ti/tj * ( (tH + s4)/tj + (ti - uj)/sH );
 551   } else {
 552     fac1 = -ti/sH + 2.0 * ( uH*tH - s4*s3 )/sH/uj;
 553     fac2 = ui/uj * ( (uH + s4)/uj + (ui - tj)/sH );
 554   }
 555
 556   // Compute matrix element weight
 557   double weight = 0.0;
 558
 559   // Average over separate helicity contributions
 560   // (for qbar g : ha -> -ha )
 561   // LL (ha = -1, hb = +1) (divided by 4 for average)
 562   weight += fac2 * norm(LsqqX) / 2.0;
 563   // RR (ha =  1, hb = -1) (divided by 4 for average)
 564   weight += fac2 * norm(RsqqX) / 2.0;
 565   // RL (ha =  1, hb =  1) (divided by 4 for average)
 566   weight += fac2 * norm(RsqqX) / 2.0 + fac1 * norm(RsqqX);
 567   // LR (ha = -1, hb = -1) (divided by 4 for average)
 568   weight += fac2 * norm(LsqqX) / 2.0 + fac1 * norm(LsqqX);
 569
 570   double sigma = sigma0 * weight;
 571   if (abs(idq) < 9) sigma /= 3.;
 572
 573   // Answer.
 574   return sigma;
 575
 576 }
 577
 578 //--------------------------------------------------------------------------
 579
 580 // Select identity, colour and anticolour.
 581
 582 void Sigma2qg2chi0squark::setIdColAcol() {
 583
 584   // Set flavours.
 585   setId( id1, id2, id3, (id1*id2 > 0 ? abs(id4) : -abs(id4)));
 586
 587   // Colour flow topology. Swap if first is gluon, or when antiquark.
 588   if (id1 != 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
 589   else setColAcol( 1, 2, 2, 0, 0, 0, 1, 0);
 590   if (id1*id2 < 0) swapColAcol();
 591
 592 }
 593
 594 //==========================================================================
 595
 596 // Sigma2qg2charsquark
 597 // Cross section for gaugino-squark production: chargino-squark
 598
 599 //--------------------------------------------------------------------------
 600
 601 // Initialize process.
 602
 603 void Sigma2qg2charsquark::initProc() {
 604
 605   //Typecast to the correct couplings
 606   coupSUSYPtr = (CoupSUSY*) couplingsPtr;
 607
 608   // Construct name of process.
 609   if (id4 % 2 == 0) {
 610     nameSave = "q g -> " + particleDataPtr->name(id3) + " "
 611       + particleDataPtr->name(id4) + " + c.c. (q=d,s,b)";
 612   }
 613   else {
 614     nameSave = "q g -> " + particleDataPtr->name(id3) + " "
 615       + particleDataPtr->name(id4) + " + c.c. (q=u,c)";
 616   }
 617
 618   // Secondary open width fraction.
 619   openFracPair = particleDataPtr->resOpenFrac(id3Sav, id4Sav);
 620
 621 }
 622
 623 //--------------------------------------------------------------------------
 624
 625 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
 626
 627 double Sigma2qg2charsquark::sigmaHat() {
 628
 629   // Antiquark -> antisquark
 630   int idq = (id1 == 21) ? id2 : id1;
 631   if (idq > 0) {
 632     id3 = id3Sav;
 633     id4 = id4Sav;
 634   } else {
 635     id3 = -id3Sav;
 636     id4 = -id4Sav;
 637   }
 638
 639   // Only accept u(bar) -> ~d(bar) and d(bar) -> ~u(bar)
 640   if (particleDataPtr->chargeType(idq) == particleDataPtr->chargeType(id4))
 641     return 0.0;
 642
 643   // Generation index
 644   int iGq = (abs(idq)+1)/2;
 645
 646   // Couplings
 647   complex LsqqX, RsqqX;
 648   if (idq % 2 == 0) {
 649     LsqqX = coupSUSYPtr->LsduX[id4sq][iGq][id3chi];
 650     RsqqX = coupSUSYPtr->RsduX[id4sq][iGq][id3chi];
 651   }
 652   else {
 653     LsqqX = coupSUSYPtr->LsduX[id4sq][iGq][id3chi];
 654     RsqqX = coupSUSYPtr->RsduX[id4sq][iGq][id3chi];
 655   }
 656
 657   // Prefactors : swap u and t if gq instead of qg
 658   double fac1, fac2;
 659   if (idq == id1) {
 660     fac1 = -ui/sH + 2.0 * ( uH*tH - s4*s3 )/sH/tj;
 661     fac2 = ti/tj * ( (tH + s4)/tj + (ti - uj)/sH );
 662   } else {
 663     fac1 = -ti/sH + 2.0 * ( uH*tH - s4*s3 )/sH/uj;
 664     fac2 = ui/uj * ( (uH + s4)/uj + (ui - tj)/sH );
 665   }
 666
 667   // Compute matrix element weight
 668   double weight = 0.0;
 669
 670   // Average over separate helicity contributions
 671   // (a, b refers to qg configuration)
 672   // LL (ha = -1, hb = +1) (divided by 4 for average)
 673   weight += fac2 * norm(LsqqX) / 2.0;
 674   // RR (ha =  1, hb = -1) (divided by 4 for average)
 675   weight += fac2 * norm(RsqqX) / 2.0;
 676   // RL (ha =  1, hb =  1) (divided by 4 for average)
 677   weight += fac2 * norm(RsqqX) / 2.0 + fac1 * norm(RsqqX);
 678   // LR (ha = -1, hb = -1) (divided by 4 for average)
 679   weight += fac2 * norm(LsqqX) / 2.0 + fac1 * norm(LsqqX);
 680
 681   double sigma = sigma0 * weight;
 682   if (abs(idq) < 9) sigma /= 3.;
 683
 684   // Answer.
 685   return sigma * openFracPair;
 686
 687 }
 688
 689 //--------------------------------------------------------------------------
 690
 691 // Select identity, colour and anticolour.
 692
 693 void Sigma2qg2charsquark::setIdColAcol() {
 694
 695   // Set flavours.
 696   if (id1 > 0 && id2 > 0) {
 697     setId( id1, id2, id3Sav, id4Sav);
 698   } else {
 699     setId( id1, id2,-id3Sav,-id4Sav);
 700   }
 701
 702   // Colour flow topology. Swap if first is gluon, or when antiquark.
 703   if (id1 != 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
 704   else setColAcol( 1, 2, 2, 0, 0, 0, 1, 0);
 705   if (id1 < 0 || id2 < 0) swapColAcol();
 706
 707 }
 708
 709 //==========================================================================
 710
 711 // Sigma2qq2squarksquark
 712 // Cross section for squark-squark production
 713
 714 //--------------------------------------------------------------------------
 715
 716 // Initialize process.
 717
 718 void Sigma2qq2squarksquark::initProc() {
 719
 720   //Typecast to the correct couplings
 721   coupSUSYPtr = (CoupSUSY*) couplingsPtr;
 722
 723   // Extract mass-ordering indices
 724   iGen3 = 3*(abs(id3Sav)/2000000) + (abs(id3Sav)%10+1)/2;
 725   iGen4 = 3*(abs(id4Sav)/2000000) + (abs(id4Sav)%10+1)/2;
 726
 727   // Is this a ~u_i ~d_j fial state or ~d_i ~d_j, ~u_i ~u_j
 728   if (abs(id3Sav) % 2 == abs(id4Sav) % 2) isUD = false;
 729   else isUD = true;
 730
 731   // Derive name
 732   nameSave = "q q' -> "+particleDataPtr->name(abs(id3Sav))+" "
 733     +particleDataPtr->name(abs(id4Sav))+" + c.c.";
 734
 735   // Count 5 neutralinos in NMSSM
 736   nNeut = (coupSUSYPtr->isNMSSM ? 5 : 4);
 737
 738   // Store mass squares of all possible internal propagator lines
 739   m2Glu     = pow2(particleDataPtr->m0(1000021));
 740   m2Neut.resize(nNeut+1);
 741   for (int iNeut=1;iNeut<=nNeut;iNeut++) {
 742     m2Neut[iNeut] = pow2(particleDataPtr->m0(coupSUSYPtr->idNeut(iNeut)));
 743   }
 744   m2Char.resize(3);
 745   m2Char[1] = pow2(particleDataPtr->m0(coupSUSYPtr->idChar(1)));
 746   m2Char[2] = pow2(particleDataPtr->m0(coupSUSYPtr->idChar(2)));
 747
 748   // Set sizes of some arrays to be used below
 749   tNeut.resize(nNeut+1);
 750   uNeut.resize(nNeut+1);
 751   tChar.resize(3);
 752   uChar.resize(3);
 753
 754   // Secondary open width fraction.
 755   openFracPair = particleDataPtr->resOpenFrac(id3Sav, id4Sav);
 756 }
 757
 758 //--------------------------------------------------------------------------
 759
 760 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
 761
 762 void Sigma2qq2squarksquark::sigmaKin() {
 763
 764   // Weak mixing
 765   double xW = coupSUSYPtr->sin2W;
 766
 767   // pi/sH2
 768   double comFacHat = M_PI/sH2 * openFracPair;
 769
 770   // Channel-dependent but flavor-independent pre-factors
 771   sigmaNeut     = comFacHat * pow2(alpEM) / pow2(xW) / pow2(1-xW);
 772   sigmaGlu      = comFacHat * 2.0 * pow2(alpS) / 9.0;
 773   if (isUD) {
 774     sigmaChar     = comFacHat * pow2(alpEM) / 4.0 / pow2(xW);
 775     sigmaCharNeut = comFacHat * pow2(alpEM) / 3.0 / pow2(xW) / (1-xW);
 776     sigmaCharGlu  = comFacHat * 4.0 * alpEM * alpS / 9.0 / xW;
 777     sigmaNeutGlu  = 0.0;
 778   } else {
 779     sigmaChar     = 0.0;
 780     sigmaCharNeut = 0.0;
 781     sigmaCharGlu  = 0.0;
 782     sigmaNeutGlu  = comFacHat * 8.0 * alpEM * alpS / 9.0 / xW/(1-xW);
 783   }
 784
 785 }
 786
 787 //--------------------------------------------------------------------------
 788
 789 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
 790
 791 double Sigma2qq2squarksquark::sigmaHat() {
 792
 793   // In-pair must be same-sign
 794   if (id1 * id2 < 0) return 0.0;
 795
 796   // Check correct charge sum
 797   if (isUD && abs(id1) %2 == abs(id2) % 2) return 0.0;
 798   if (!isUD && abs(id1) % 2 != abs(id2) % 2) return 0.0;
 799   if (!isUD && abs(id1) % 2 != abs(id3Sav) % 2) return 0.0;
 800
 801   // Coded sigma is for ud -> ~q~q'. Swap t and u for du -> ~q~q'.
 802   swapTU = (isUD and abs(id1) % 2 == 0);
 803   int    idIn1A = (swapTU) ? abs(id2) : abs(id1);
 804   int    idIn2A = (swapTU) ? abs(id1) : abs(id2);
 805
 806   // Auxiliary factors for use below
 807   tGlu     = tH - m2Glu;
 808   uGlu     = uH - m2Glu;
 809   for (int i=1; i<= nNeut; i++) {
 810     tNeut[i] = tH - m2Neut[i];
 811     uNeut[i] = uH - m2Neut[i];
 812     if (isUD && i <= 2) {
 813       tChar[i] = tH - m2Char[i];
 814       uChar[i] = uH - m2Char[i];
 815     }
 816   }
 817
 818   // Generation indices of incoming particles
 819   int iGen1 = (abs(idIn1A)+1)/2;
 820   int iGen2 = (abs(idIn2A)+1)/2;
 821
 822   // Initial values for pieces used for color-flow selection below
 823   sumCt = 0.0;
 824   sumCu = 0.0;
 825   sumNt = 0.0;
 826   sumNu = 0.0;
 827   sumGt = 0.0;
 828   sumGu = 0.0;
 829   sumInterference = 0.0;
 830
 831   // Common factor for LR and RL contributions
 832   double facTU =  uH*tH-s3*s4;
 833
 834   // Case A) Opposite-isospin: qq' -> ~d~u
 835   if ( isUD ) {
 836
 837     // t-channel Charginos
 838     for (int k=1;k<=2;k++) {
 839
 840       // Skip if only including gluinos
 841       if (settingsPtr->flag("SUSY:qq2squarksquark:onlyQCD")) continue;
 842
 843       for (int l=1;l<=2;l++) {
 844
 845         // kl-dependent factor for LL and RR contributions
 846         double facMS = sH*sqrt(m2Char[k]*m2Char[l]);
 847
 848         // Note: Ckl defined as in [Boz07] with sigmaChar factored out
 849         // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
 850         complex Ckl[3][3];
 851         Ckl[1][1] = facMS
 852           * coupSUSYPtr->LsudX[iGen4][iGen2][k]
 853           * conj(coupSUSYPtr->LsduX[iGen3][iGen1][k])
 854           * conj(coupSUSYPtr->LsudX[iGen4][iGen2][l])
 855           * coupSUSYPtr->LsduX[iGen3][iGen1][l];
 856         Ckl[1][2] = facTU
 857           * coupSUSYPtr->RsudX[iGen4][iGen2][k]
 858           * conj(coupSUSYPtr->LsduX[iGen3][iGen1][k])
 859           * conj(coupSUSYPtr->RsudX[iGen4][iGen2][l])
 860           * coupSUSYPtr->LsduX[iGen3][iGen1][l];
 861         Ckl[2][1] = facTU
 862           * coupSUSYPtr->LsudX[iGen4][iGen2][k]
 863           * conj(coupSUSYPtr->RsduX[iGen3][iGen1][k])
 864           * conj(coupSUSYPtr->LsudX[iGen4][iGen2][l])
 865           * coupSUSYPtr->RsduX[iGen3][iGen1][l];
 866         Ckl[2][2] = facMS
 867           * coupSUSYPtr->RsudX[iGen4][iGen2][k]
 868           * conj(coupSUSYPtr->RsduX[iGen3][iGen1][k])
 869           * conj(coupSUSYPtr->RsudX[iGen4][iGen2][l])
 870           * coupSUSYPtr->RsduX[iGen3][iGen1][l];
 871
 872         // Add to sum of t-channel charginos
 873         sumCt += sigmaChar * real(Ckl[1][1] + Ckl[1][2] + Ckl[2][1] + Ckl[2][2])
 874           / tChar[k] / tChar[l];
 875
 876       }
 877     }
 878
 879     // u-channel Neutralinos
 880     for (int k=1;k<=nNeut;k++) {
 881
 882       // Skip if only including gluinos
 883       if (settingsPtr->flag("SUSY:qq2squarksquark:onlyQCD")) continue;
 884
 885       for (int l=1;l<=nNeut;l++) {
 886
 887         // kl-dependent factor for LL, RR contributions
 888         double facMS = sH*sqrt(m2Neut[k]*m2Neut[l]);
 889
 890         // Note: Nkl defined as in [Boz07] with sigmaNeut factored out
 891         // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
 892         complex Nkl[3][3];
 893         Nkl[1][1] = facMS
 894           * conj(coupSUSYPtr->LsuuX[iGen4][iGen1][k])
 895           * conj(coupSUSYPtr->LsddX[iGen3][iGen2][k])
 896           * coupSUSYPtr->LsuuX[iGen4][iGen1][l]
 897           * coupSUSYPtr->LsddX[iGen3][iGen2][l];
 898         Nkl[1][2] = facTU
 899           * conj(coupSUSYPtr->LsuuX[iGen4][iGen1][k])
 900           * conj(coupSUSYPtr->RsddX[iGen3][iGen2][k])
 901           * coupSUSYPtr->LsuuX[iGen4][iGen1][l]
 902           * coupSUSYPtr->RsddX[iGen3][iGen2][l];
 903         Nkl[2][1] =  facTU
 904           * conj(coupSUSYPtr->RsuuX[iGen4][iGen1][k])
 905           * conj(coupSUSYPtr->LsddX[iGen3][iGen2][k])
 906           * coupSUSYPtr->RsuuX[iGen4][iGen1][l]
 907           * coupSUSYPtr->LsddX[iGen3][iGen2][l];
 908         Nkl[2][2] =  facMS
 909           * conj(coupSUSYPtr->RsuuX[iGen4][iGen1][k])
 910           * conj(coupSUSYPtr->RsddX[iGen3][iGen2][k])
 911           * coupSUSYPtr->RsuuX[iGen4][iGen1][l]
 912           * coupSUSYPtr->RsddX[iGen3][iGen2][l];
 913
 914         // Add to sum of u-channel neutralinos
 915         sumNu += sigmaNeut / uNeut[k] / uNeut[l]
 916           * real(Nkl[1][1] + Nkl[1][2] + Nkl[2][1] + Nkl[2][2]);
 917
 918       }
 919     }
 920
 921     // u-channel gluino
 922     // Note: Gij defined as in [Boz07] with sigmaGlu factored out
 923     // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
 924     double Gij[3][3];
 925     Gij[1][1] = norm(coupSUSYPtr->LsuuG[iGen4][iGen1]
 926                 * coupSUSYPtr->LsddG[iGen3][iGen2]);
 927     Gij[1][2] = norm(coupSUSYPtr->LsuuG[iGen4][iGen1]
 928                 * coupSUSYPtr->RsddG[iGen3][iGen2]);
 929     Gij[2][1] = norm(coupSUSYPtr->RsuuG[iGen4][iGen1]
 930                 * coupSUSYPtr->LsddG[iGen3][iGen2]);
 931     Gij[2][2] = norm(coupSUSYPtr->RsuuG[iGen4][iGen1]
 932                 * coupSUSYPtr->RsddG[iGen3][iGen2]);
 933     Gij[1][1] *= sH*m2Glu;
 934     Gij[1][2] *= facTU;
 935     Gij[2][1] *= facTU;
 936     Gij[2][2] *= sH*m2Glu;
 937     // Sum over polarizations
 938     sumGu += sigmaGlu * (Gij[1][1] + Gij[1][2] + Gij[2][1] + Gij[2][2])
 939       / pow2(uGlu);
 940
 941     // chargino-neutralino interference
 942     for (int k=1;k<=2;k++) {
 943       for (int l=1;l<=nNeut;l++) {
 944
 945         // Skip if only including gluinos
 946         if (settingsPtr->flag("SUSY:qq2squarksquark:onlyQCD")) continue;
 947
 948         // Note: CNkl defined as in [Boz07] with pi/sH2 factored out
 949         // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
 950         double CNkl[3][3];
 951         CNkl[1][1] = real(coupSUSYPtr->LsudX[iGen4][iGen2][k]
 952                           * conj(coupSUSYPtr->LsduX[iGen3][iGen1][k])
 953                           * coupSUSYPtr->LsuuX[iGen4][iGen1][l]
 954                           * coupSUSYPtr->LsddX[iGen3][iGen2][l]);
 955         CNkl[1][2] = real(coupSUSYPtr->RsudX[iGen4][iGen2][k]
 956                           * conj(coupSUSYPtr->LsduX[iGen3][iGen1][k])
 957                           * coupSUSYPtr->LsuuX[iGen4][iGen1][l]
 958                           * coupSUSYPtr->RsddX[iGen3][iGen2][l]);
 959         CNkl[2][1] = real(coupSUSYPtr->LsudX[iGen4][iGen2][k]
 960                           * conj(coupSUSYPtr->RsduX[iGen3][iGen1][k])
 961                           * coupSUSYPtr->RsuuX[iGen4][iGen1][l]
 962                           * coupSUSYPtr->LsddX[iGen3][iGen2][l]);
 963         CNkl[2][2] = real(coupSUSYPtr->RsudX[iGen4][iGen2][k]
 964                           * conj(coupSUSYPtr->RsduX[iGen3][iGen1][k])
 965                           * coupSUSYPtr->RsuuX[iGen4][iGen1][l]
 966                           * coupSUSYPtr->RsddX[iGen3][iGen2][l]);
 967         CNkl[1][1] *= sH*sqrt(m2Char[k]*m2Neut[l]);
 968         CNkl[1][2] *= uH*tH-s3*s4;
 969         CNkl[2][1] *= uH*tH-s3*s4;
 970         CNkl[2][2] *= sH*sqrt(m2Char[k]*m2Neut[l]);
 971         // Sum over polarizations
 972         sumInterference += sigmaCharNeut * (CNkl[1][1] + CNkl[1][2] + CNkl[2][1]
 973                                             + CNkl[2][2]) / tChar[k] / uNeut[l];
 974       }
 975     }
 976
 977     // chargino-gluino interference
 978     for (int k=1;k<=2;k++) {
 979
 980       // Skip if only including gluinos
 981       if (settingsPtr->flag("SUSY:qq2squarksquark:onlyQCD")) continue;
 982
 983       // Note: CGk defined as in [Boz07] with sigmaCharGlu factored out
 984       // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
 985       double CGk[3][3];
 986       CGk[1][1] = real(coupSUSYPtr->LsudX[iGen4][iGen2][k]
 987                        * conj(coupSUSYPtr->LsduX[iGen3][iGen1][k])
 988                        * conj(coupSUSYPtr->LsuuG[iGen4][iGen1])
 989                        * conj(coupSUSYPtr->LsddG[iGen3][iGen2]));
 990       CGk[1][2] = real(coupSUSYPtr->RsudX[iGen4][iGen2][k]
 991                        * conj(coupSUSYPtr->LsduX[iGen3][iGen1][k])
 992                        * conj(coupSUSYPtr->LsuuG[iGen4][iGen1])
 993                        * conj(coupSUSYPtr->RsddG[iGen3][iGen2]));
 994       CGk[2][1] = real(coupSUSYPtr->LsudX[iGen4][iGen2][k]
 995                        * conj(coupSUSYPtr->RsduX[iGen3][iGen1][k])
 996                        * conj(coupSUSYPtr->RsuuG[iGen4][iGen1])
 997                        * conj(coupSUSYPtr->LsddG[iGen3][iGen2]));
 998       CGk[2][2] = real(coupSUSYPtr->RsudX[iGen4][iGen2][k]
 999                        * conj(coupSUSYPtr->RsduX[iGen3][iGen1][k])
1000                        * conj(coupSUSYPtr->RsuuG[iGen4][iGen1])
1001                        * conj(coupSUSYPtr->RsddG[iGen3][iGen2]));
1002       CGk[1][1] *= sH*sqrt(m2Glu*m2Char[k]);
1003       CGk[1][2] *= uH*tH-s3*s4;
1004       CGk[2][1] *= uH*tH-s3*s4;
1005       CGk[2][2] *= sH*sqrt(m2Glu*m2Char[k]);
1006       // Sum over polarizations
1007       sumInterference += sigmaGlu * (CGk[1][1] + CGk[1][2] + CGk[2][1]
1008         + CGk[2][2]) / uGlu / tChar[k];
1009     }
1010   }
1011
1012   // Case B) Same-isospin: qq' -> ~d~d , ~u~u
1013   else {
1014
1015     // t-channel + u-channel Neutralinos + t/u interference
1016     for (int k=1;k<=nNeut;k++) {
1017
1018       // Skip if only including gluinos
1019       if (settingsPtr->flag("SUSY:qq2squarksquark:onlyQCD")) continue;
1020
1021       for (int l=1;l<=nNeut;l++) {
1022
1023         // kl-dependent factor for LL and RR contributions
1024         double facMS = sH * particleDataPtr->m0(coupSUSYPtr->idNeut(k))
1025           * particleDataPtr->m0(coupSUSYPtr->idNeut(l));
1026
1027         // Note: Nxkl defined as in [Boz07] with sigmaNeut factored out
1028         // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
1029         complex NTkl[3][3], NUkl[3][3], NTUkl[3][3];
1030         NTkl[1][1] = facMS
1031           * conj(coupSUSYPtr->getLsqqX(iGen4,idIn2A,k))
1032           * conj(coupSUSYPtr->getLsqqX(iGen3,idIn1A,k))
1033           * coupSUSYPtr->getLsqqX(iGen4,idIn2A,l)
1034           * coupSUSYPtr->getLsqqX(iGen3,idIn1A,l);
1035         NTkl[1][2] = facTU
1036           * conj(coupSUSYPtr->getRsqqX(iGen4,idIn2A,k))
1037           * conj(coupSUSYPtr->getLsqqX(iGen3,idIn1A,k))
1038           * coupSUSYPtr->getRsqqX(iGen4,idIn2A,l)
1039           * coupSUSYPtr->getLsqqX(iGen3,idIn1A,l);
1040         NTkl[2][1] = facTU
1041           * conj(coupSUSYPtr->getLsqqX(iGen4,idIn2A,k))
1042           * conj(coupSUSYPtr->getRsqqX(iGen3,idIn1A,k))
1043           * coupSUSYPtr->getLsqqX(iGen4,idIn2A,l)
1044           * coupSUSYPtr->getRsqqX(iGen3,idIn1A,l);
1045         NTkl[2][2] = facMS
1046           * conj(coupSUSYPtr->getRsqqX(iGen4,idIn2A,k))
1047           * conj(coupSUSYPtr->getRsqqX(iGen3,idIn1A,k))
1048           * coupSUSYPtr->getRsqqX(iGen4,idIn2A,l)
1049           * coupSUSYPtr->getRsqqX(iGen3,idIn1A,l);
1050         NUkl[1][1] = facMS
1051           * conj(coupSUSYPtr->getLsqqX(iGen3,idIn2A,k))
1052           * conj(coupSUSYPtr->getLsqqX(iGen4,idIn1A,k))
1053           * coupSUSYPtr->getLsqqX(iGen3,idIn2A,l)
1054           * coupSUSYPtr->getLsqqX(iGen4,idIn1A,l);
1055         NUkl[1][2] = facTU
1056           * conj(coupSUSYPtr->getRsqqX(iGen3,idIn2A,k))
1057           * conj(coupSUSYPtr->getLsqqX(iGen4,idIn1A,k))
1058           * coupSUSYPtr->getRsqqX(iGen3,idIn2A,l)
1059           * coupSUSYPtr->getLsqqX(iGen4,idIn1A,l);
1060         NUkl[2][1] = facTU
1061           * conj(coupSUSYPtr->getLsqqX(iGen3,idIn2A,k))
1062           * conj(coupSUSYPtr->getRsqqX(iGen4,idIn1A,k))
1063           * coupSUSYPtr->getLsqqX(iGen3,idIn2A,l)
1064           * coupSUSYPtr->getRsqqX(iGen4,idIn1A,l);
1065         NUkl[2][2] = facMS
1066           * conj(coupSUSYPtr->getRsqqX(iGen3,idIn2A,k))
1067           * conj(coupSUSYPtr->getRsqqX(iGen4,idIn1A,k))
1068           * coupSUSYPtr->getRsqqX(iGen3,idIn2A,l)
1069           * coupSUSYPtr->getRsqqX(iGen4,idIn1A,l);
1070         NTUkl[1][1] = facMS
1071           * real( conj(coupSUSYPtr->getLsqqX(iGen4,idIn2A,k))
1072                   * conj(coupSUSYPtr->getLsqqX(iGen3,idIn1A,k))
1073                   * coupSUSYPtr->getLsqqX(iGen3,idIn2A,l)
1074                   * coupSUSYPtr->getLsqqX(iGen4,idIn1A,l) );
1075         NTUkl[1][2] = facTU
1076           * real( conj(coupSUSYPtr->getRsqqX(iGen4,idIn2A,k))
1077                   * conj(coupSUSYPtr->getLsqqX(iGen3,idIn1A,k))
1078                   * coupSUSYPtr->getRsqqX(iGen3,idIn2A,l)
1079                   * coupSUSYPtr->getLsqqX(iGen4,idIn1A,l) );
1080         NTUkl[2][1] = facTU
1081           * real( conj(coupSUSYPtr->getLsqqX(iGen4,idIn2A,k))
1082                   * conj(coupSUSYPtr->getRsqqX(iGen3,idIn1A,k))
1083                   * coupSUSYPtr->getLsqqX(iGen3,idIn2A,l)
1084                   * coupSUSYPtr->getRsqqX(iGen4,idIn1A,l) );
1085         NTUkl[2][2] = facMS
1086           * real( conj(coupSUSYPtr->getRsqqX(iGen4,idIn2A,k))
1087                   * conj(coupSUSYPtr->getRsqqX(iGen3,idIn1A,k))
1088                   * coupSUSYPtr->getRsqqX(iGen3,idIn2A,l)
1089                   * coupSUSYPtr->getRsqqX(iGen4,idIn1A,l) );
1090
1091         // Add to sums
1092         sumNt += sigmaNeut / tNeut[k] / tNeut[l]
1093           * real(NTkl[1][1] + NTkl[1][2] + NTkl[2][1] + NTkl[2][2]);
1094         sumNu += sigmaNeut / uNeut[k] / uNeut[l]
1095           * real(NUkl[1][1] + NUkl[1][2] + NUkl[2][1] + NUkl[2][2]);
1096         sumInterference += 2.0 / 3.0 * sigmaNeut
1097           * real(NTUkl[1][1] + NTUkl[1][2] + NTUkl[2][1] + NTUkl[2][2])
1098           / tNeut[k] / uNeut[l];
1099       }
1100
1101       // Neutralino / Gluino interference
1102
1103       // k-dependent factor for LL and RR contributions
1104       double facMS = sH * particleDataPtr->m0(coupSUSYPtr->idNeut(k))
1105         * particleDataPtr->m0(1000021);
1106
1107       // Note: Nxkl defined as in [Boz07] with sigmaNeutGlu factored out
1108       // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
1109       complex NGA[3][3], NGB[3][3];
1110       NGA[1][1] = facMS
1111         * real( conj(coupSUSYPtr->getLsqqX(iGen4,idIn2A,k))
1112                 * conj(coupSUSYPtr->getLsqqX(iGen3,idIn1A,k))
1113                 * conj(coupSUSYPtr->getLsqqG(iGen3,idIn2A))
1114                 * conj(coupSUSYPtr->getLsqqG(iGen4,idIn1A)) );
1115       NGA[1][2] = facTU
1116         * real( conj(coupSUSYPtr->getRsqqX(iGen4,idIn2A,k))
1117                 * conj(coupSUSYPtr->getLsqqX(iGen3,idIn1A,k))
1118                 * conj(coupSUSYPtr->getLsqqG(iGen3,idIn2A))
1119                 * conj(coupSUSYPtr->getRsqqG(iGen4,idIn1A)) );
1120       NGA[2][1] = facTU
1121         * real( conj(coupSUSYPtr->getLsqqX(iGen4,idIn2A,k))
1122                 * conj(coupSUSYPtr->getRsqqX(iGen3,idIn1A,k))
1123                 * conj(coupSUSYPtr->getRsqqG(iGen3,idIn2A))
1124                 * conj(coupSUSYPtr->getLsqqG(iGen4,idIn1A)) );
1125       NGA[2][2] = facMS
1126         * real( conj(coupSUSYPtr->getRsqqX(iGen4,idIn2A,k))
1127                 * conj(coupSUSYPtr->getRsqqX(iGen3,idIn1A,k))
1128                 * conj(coupSUSYPtr->getRsqqG(iGen3,idIn2A))
1129                 * conj(coupSUSYPtr->getRsqqG(iGen4,idIn1A)) );
1130       NGB[1][1] = facMS
1131         * real( conj(coupSUSYPtr->getLsqqX(iGen3,idIn2A,k))
1132                 * conj(coupSUSYPtr->getLsqqX(iGen4,idIn1A,k))
1133                 * conj(coupSUSYPtr->getLsqqG(iGen4,idIn2A))
1134                 * conj(coupSUSYPtr->getLsqqG(iGen3,idIn1A)) );
1135       NGB[1][2] = facMS
1136         * real( conj(coupSUSYPtr->getRsqqX(iGen3,idIn2A,k))
1137                 * conj(coupSUSYPtr->getLsqqX(iGen4,idIn1A,k))
1138                 * conj(coupSUSYPtr->getRsqqG(iGen4,idIn2A))
1139                 * conj(coupSUSYPtr->getLsqqG(iGen3,idIn1A)) );
1140       NGB[2][1] = facMS
1141         * real( conj(coupSUSYPtr->getLsqqX(iGen3,idIn2A,k))
1142                 * conj(coupSUSYPtr->getRsqqX(iGen4,idIn1A,k))
1143                 * conj(coupSUSYPtr->getLsqqG(iGen4,idIn2A))
1144                 * conj(coupSUSYPtr->getRsqqG(iGen3,idIn1A)) );
1145       NGB[2][2] = facMS
1146         * real( conj(coupSUSYPtr->getRsqqX(iGen3,idIn2A,k))
1147                 * conj(coupSUSYPtr->getRsqqX(iGen4,idIn1A,k))
1148                 * conj(coupSUSYPtr->getRsqqG(iGen4,idIn2A))
1149                 * conj(coupSUSYPtr->getRsqqG(iGen3,idIn1A)) );
1150
1151       // Add to sums
1152       sumInterference += sigmaNeutGlu *
1153         ( real(NGA[1][1] + NGA[1][2] + NGA[2][1] + NGA[2][2]) / tNeut[k] / uGlu
1154         + real(NGB[1][1] + NGB[1][2] + NGB[2][1] + NGB[2][2]) / uNeut[k] / tGlu );
1155     }
1156
1157     // t-channel + u-channel Gluinos + t/u interference
1158
1159     // factor for LL and RR contributions
1160     double facMS = sH * m2Glu;
1161
1162     // Note: GT, GU defined as in [Boz07] with sigmaGlu factored out
1163     // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
1164     complex GT[3][3], GU[3][3], GTU[3][3];
1165     GT[1][1] = facMS
1166       * norm(coupSUSYPtr->getLsqqG(iGen4,idIn2A)
1167       * coupSUSYPtr->getLsqqG(iGen3,idIn1A));
1168     GT[1][2] = facTU
1169       * norm(coupSUSYPtr->getRsqqG(iGen4,idIn2A)
1170       * coupSUSYPtr->getLsqqG(iGen3,idIn1A));
1171     GT[2][1] = facTU
1172       * norm(coupSUSYPtr->getLsqqG(iGen4,idIn2A)
1173       * coupSUSYPtr->getRsqqG(iGen3,idIn1A));
1174     GT[2][2] = facMS
1175       * norm(coupSUSYPtr->getRsqqG(iGen4,idIn2A)
1176       * coupSUSYPtr->getRsqqG(iGen3,idIn1A));
1177     GU[1][1] = facMS
1178       * norm(coupSUSYPtr->getLsqqG(iGen3,idIn2A)
1179       * coupSUSYPtr->getLsqqG(iGen4,idIn1A));
1180     GU[1][2] = facTU
1181       * norm(coupSUSYPtr->getLsqqG(iGen3,idIn2A)
1182       * coupSUSYPtr->getRsqqG(iGen4,idIn1A));
1183     GU[2][1] = facTU
1184       * norm(coupSUSYPtr->getRsqqG(iGen3,idIn2A)
1185       * coupSUSYPtr->getLsqqG(iGen4,idIn1A));
1186     GU[2][2] = facMS
1187       * norm(coupSUSYPtr->getRsqqG(iGen3,idIn2A)
1188       * coupSUSYPtr->getRsqqG(iGen4,idIn1A));
1189
1190     GTU[1][1] = facMS
1191       * real(coupSUSYPtr->getLsqqG(iGen3,idIn1A)
1192              * coupSUSYPtr->getLsqqG(iGen4,idIn2A)
1193              * conj(coupSUSYPtr->getLsqqG(iGen3,idIn2A))
1194              * conj(coupSUSYPtr->getLsqqG(iGen4,idIn1A)) );
1195
1196     GTU[1][2] = facTU
1197       * real(coupSUSYPtr->getLsqqG(iGen3,idIn1A)
1198              * coupSUSYPtr->getRsqqG(iGen4,idIn2A)
1199              * conj(coupSUSYPtr->getRsqqG(iGen3,idIn2A))
1200              * conj(coupSUSYPtr->getLsqqG(iGen4,idIn1A)) );
1201
1202     GTU[2][1] = facTU
1203       * real(coupSUSYPtr->getRsqqG(iGen3,idIn1A)
1204              * coupSUSYPtr->getLsqqG(iGen4,idIn2A)
1205              * conj(coupSUSYPtr->getLsqqG(iGen3,idIn2A))
1206              * conj(coupSUSYPtr->getRsqqG(iGen4,idIn1A)) );
1207
1208     GTU[2][2] = facMS
1209       * real(coupSUSYPtr->getRsqqG(iGen3,idIn1A)
1210              * coupSUSYPtr->getRsqqG(iGen4,idIn2A)
1211              * conj(coupSUSYPtr->getRsqqG(iGen3,idIn2A))
1212              * conj(coupSUSYPtr->getRsqqG(iGen4,idIn1A)) );
1213
1214     // Add to sums
1215     sumGt += sigmaGlu * real(GT[1][1] + GT[1][2] + GT[2][1] + GT[2][2])
1216       / pow2(tGlu) ;
1217     sumGu += sigmaGlu * real(GU[1][1] + GU[1][2] + GU[2][1] + GU[2][2])
1218       / pow2(uGlu) ;
1219     sumInterference += - 2.0 / 3.0 * sigmaGlu
1220       * real(GTU[1][1] + GTU[1][2] + GTU[2][1] + GTU[2][2])
1221       / tGlu / uGlu;
1222
1223   }
1224
1225   // Cross section
1226   double sigma = sumNt + sumNu + sumCt + sumCu + sumGt + sumGu
1227     + sumInterference;
1228
1229   // Identical particles?
1230   if (id3Sav == id4Sav) sigma /= 2.0;
1231
1232   // Return answer.
1233   return sigma;
1234
1235 }
1236
1237 //--------------------------------------------------------------------------
1238
1239 // Select identity, colour and anticolour.
1240
1241 void Sigma2qq2squarksquark::setIdColAcol() {
1242
1243   // Set flavours.
1244   if (id1 > 0 && id2 > 0) {
1245     setId( id1, id2, id3Sav, id4Sav);
1246   } else {
1247     // 1,2 -> -3,-4
1248     setId( id1, id2,-id3Sav,-id4Sav);
1249   }
1250
1251   // Coded sigma is for ud -> ~q~q'. Swap t and u for du -> ~q~q'.
1252   swapTU = (isUD and abs(id1) % 2 == 0);
1253
1254   // Select colour flow topology
1255   // A: t-channel neutralino, t-channel chargino, or u-channel gluino
1256   double fracA = sumNt + sumCt + sumGu
1257     / (sumNt + sumNu + sumCt + sumCu + sumGt + sumGu);
1258   if (swapTU) fracA = 1.0 - fracA;
1259   setColAcol( 1, 0, 2, 0, 1, 0, 2, 0);
1260   // B: t-channel gluino or u-channel neutralino
1261   if (rndmPtr->flat() > fracA) setColAcol( 1, 0, 2, 0, 2, 0, 1, 0);
1262
1263   // Switch to anti-colors if antiquarks
1264   if (id1 < 0 || id2 < 0) swapColAcol();
1265
1266 }
1267
1268 //==========================================================================
1269
1270 // Sigma2qqbar2squarkantisquark
1271 // Cross section for qqbar-initiated squark-antisquark production
1272
1273 //--------------------------------------------------------------------------
1274
1275 // Initialize process.
1276
1277 void Sigma2qqbar2squarkantisquark::initProc() {
1278
1279   //Typecast to the correct couplings
1280   coupSUSYPtr = (CoupSUSY*) couplingsPtr;
1281
1282   // Extract isospin and mass-ordering indices
1283   iGen3 = 3*(abs(id3Sav)/2000000) + (abs(id3Sav)%10+1)/2;
1284   iGen4 = 3*(abs(id4Sav)/2000000) + (abs(id4Sav)%10+1)/2;
1285
1286   // Is this a ~u_i ~d*_j, ~d_i ~u*_j final state or ~d_i ~d*_j, ~u_i ~u*_j
1287   if (abs(id3Sav) % 2 == abs(id4Sav) % 2) isUD = false;
1288   else isUD = true;
1289
1290   // Derive name
1291   nameSave = "q qbar' -> "+particleDataPtr->name(abs(id3Sav))+" "+
1292     particleDataPtr->name(-abs(id4Sav));
1293   if (isUD && abs(id3Sav) != abs(id4Sav)) nameSave +=" + c.c.";
1294
1295   // Count 5 neutralinos in NMSSM
1296   nNeut = (coupSUSYPtr->isNMSSM ? 5 : 4);
1297
1298   // Store mass squares of all possible internal propagator lines
1299   m2Glu     = pow2(particleDataPtr->m0(1000021));
1300   m2Neut.resize(nNeut+1);
1301   for (int iNeut=1;iNeut<=nNeut;iNeut++)
1302     m2Neut[iNeut] = pow2(particleDataPtr->m0(coupSUSYPtr->idNeut(iNeut)));
1303
1304   // Set sizes of some arrays to be used below
1305   tNeut.resize(nNeut+1);
1306   uNeut.resize(nNeut+1);
1307
1308   // Shorthand for Weak mixing
1309   xW = coupSUSYPtr->sin2W;
1310
1311   // Secondary open width fraction.
1312   openFracPair = particleDataPtr->resOpenFrac(id3Sav, id4Sav);
1313
1314 }
1315
1316 //--------------------------------------------------------------------------
1317
1318 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
1319
1320 void Sigma2qqbar2squarkantisquark::sigmaKin() {
1321
1322   // Z/W propagator
1323   if (! isUD) {
1324     double sV= sH - pow2(coupSUSYPtr->mZpole);
1325     double d = pow2(sV) + pow2(coupSUSYPtr->mZpole * coupSUSYPtr->wZpole);
1326     propZW   = complex( sV / d, coupSUSYPtr->mZpole * coupSUSYPtr->wZpole / d);
1327   } else {
1328     double sV= sH - pow2(coupSUSYPtr->mWpole);
1329     double d = pow2(sV) + pow2(coupSUSYPtr->mWpole * coupSUSYPtr->wWpole);
1330     propZW   = complex( sV / d, coupSUSYPtr->mWpole * coupSUSYPtr->wWpole / d);
1331   }
1332
1333   // Flavor-independent pre-factors
1334   double comFacHat = M_PI/sH2 * openFracPair;
1335
1336   sigmaEW       = comFacHat * pow2(alpEM);
1337   sigmaGlu      = comFacHat * 2.0 * pow2(alpS) / 9.0;
1338   sigmaEWG      = comFacHat * 8.0 * alpEM * alpS / 9.0;
1339
1340 }
1341
1342 //--------------------------------------------------------------------------
1343
1344 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1345
1346 double Sigma2qqbar2squarkantisquark::sigmaHat() {
1347
1348   // In-pair must be opposite-sign
1349   if (id1 * id2 > 0) return 0.0;
1350
1351   // Check correct charge sum
1352   if (isUD && abs(id1) %2 == abs(id2) % 2) return 0.0;
1353   if (!isUD && abs(id1) % 2 != abs(id2) % 2) return 0.0;
1354
1355   // Check if using QCD diagrams only
1356   bool onlyQCD = settingsPtr->flag("SUSY:qqbar2squarkantisquark:onlyQCD");
1357
1358   // Coded UD sigma is for udbar -> ~u~d'*. Swap t<->u for dbaru -> ~u~d'*.
1359   swapTU = (isUD and abs(id1) % 2 != 0);
1360
1361   // Coded QQ sigma is for qqbar -> ~q~q*. Swap t<->u for qbarq -> ~q~q*.
1362   if (!isUD && id1 < 0) swapTU = true;
1363
1364   // Generation indices of incoming particles
1365   int idIn1A = (swapTU) ? abs(id2) : abs(id1);
1366   int idIn2A = (swapTU) ? abs(id1) : abs(id2);
1367   int iGen1  = (idIn1A+1)/2;
1368   int iGen2  = (idIn2A+1)/2;
1369
1370   // Auxiliary factors for use below
1371   tGlu     = tH - m2Glu;
1372   uGlu     = uH - m2Glu;
1373   for (int i=1; i<= nNeut; i++) {
1374     tNeut[i] = tH - m2Neut[i];
1375     uNeut[i] = uH - m2Neut[i];
1376   }
1377
1378   // Initial values for pieces used for color-flow selection below
1379   sumColS   = 0.0;
1380   sumColT   = 0.0;
1381   sumInterference = 0.0;
1382
1383   // Common factor for LR and RL contributions
1384   double facTU =  uH*tH-s3*s4;
1385
1386   // Case A) Opposite-isospin: udbar -> ~u~d*
1387   if ( isUD ) {
1388
1389     // s-channel W contribution (only contributes to LL helicities)
1390     if ( !onlyQCD ) {
1391       sumColS += sigmaEW / 16.0 / pow2(xW) / pow2(1.0-xW)
1392         * norm(conj(coupSUSYPtr->LudW[iGen1][iGen2])
1393                * coupSUSYPtr->LsusdW[iGen3][iGen4]) * facTU
1394         * norm(propZW);
1395     }
1396
1397     // t-channel gluino contributions
1398     double GT[3][3];
1399     double facLR = m2Glu * sH;
1400     // LL, LR, RL, RR
1401     GT[1][1] = facTU * norm(conj(coupSUSYPtr->LsddG[iGen4][iGen2])
1402                             *coupSUSYPtr->LsuuG[iGen3][iGen1]);
1403     GT[1][2] = facLR * norm(conj(coupSUSYPtr->RsddG[iGen4][iGen2])
1404                             *coupSUSYPtr->LsuuG[iGen3][iGen1]);
1405     GT[2][1] = facLR * norm(conj(coupSUSYPtr->LsddG[iGen4][iGen2])
1406                             *coupSUSYPtr->RsuuG[iGen3][iGen1]);
1407     GT[2][2] = facTU * norm(conj(coupSUSYPtr->RsddG[iGen4][iGen2])
1408                             *coupSUSYPtr->RsuuG[iGen3][iGen1]);
1409     // leading color flow for t-channel gluino is annihilation-like
1410     sumColS += sigmaGlu / pow2(tGlu)
1411       * (GT[1][1] + GT[1][2] + GT[2][1] + GT[2][2]);
1412
1413     // W-Gluino interference (only contributes to LL helicities)
1414     if ( !onlyQCD ) {
1415       sumColS += sigmaEWG / 4.0 / xW / (1-xW)
1416         * real(conj(coupSUSYPtr->LsuuG[iGen3][iGen1])
1417                * coupSUSYPtr->LsddG[iGen4][iGen2]
1418                * conj(coupSUSYPtr->LudW[iGen1][iGen2])
1419                * coupSUSYPtr->LsusdW[iGen3][iGen4]) * facTU
1420         / tGlu * sqrt(norm(propZW));
1421     }
1422
1423     // t-channel neutralinos
1424     // NOT YET IMPLEMENTED !
1425
1426   }
1427
1428   // Case B) Same-isospin: qqbar -> ~d~d* , ~u~u*
1429   else {
1430
1431     double eQ  = (idIn1A % 2 == 0) ? 2./3. : 1./3. ;
1432     double eSq = (abs(id3Sav) % 2 == 0) ? 2./3. : 1./3. ;
1433
1434     // s-channel gluon (strictly flavor-diagonal)
1435     if (abs(id3Sav) == abs(id4Sav) && abs(id1) == abs(id2)) {
1436       // Factor 2 since contributes to both ha != hb helicities
1437       sumColT += 2. * sigmaGlu * facTU / pow2(sH);
1438     }
1439
1440     // t-channel gluino (only for in-isospin = out-isospin).
1441     if (eQ == eSq) {
1442       // Sum over helicities.
1443       double GT[3][3];
1444       double facLR = sH * m2Glu;
1445       GT[1][1] = facTU * norm(coupSUSYPtr->getLsqqG(iGen3,idIn1A)
1446                               * conj(coupSUSYPtr->getLsqqG(iGen4,idIn2A)));
1447       GT[1][2] = facLR * norm(coupSUSYPtr->getLsqqG(iGen3,idIn1A)
1448                               * conj(coupSUSYPtr->getRsqqG(iGen4,idIn2A)));
1449       GT[2][1] = facLR * norm(coupSUSYPtr->getRsqqG(iGen3,idIn1A)
1450                               * conj(coupSUSYPtr->getLsqqG(iGen4,idIn2A)));
1451       GT[2][2] = facTU * norm(coupSUSYPtr->getRsqqG(iGen3,idIn1A)
1452                               * conj(coupSUSYPtr->getRsqqG(iGen4,idIn2A)));
1453       // Add contribution to color topology: S
1454       sumColS += sigmaGlu / pow2(tGlu)
1455         * ( GT[1][1] + GT[1][2] + GT[2][1] + GT[2][2]);
1456
1457       // gluon-gluino interference (strictly flavor-diagonal)
1458       if (abs(id3Sav) == abs(id4Sav) && abs(id1) == abs(id2)) {
1459         double GG11, GG22;
1460         GG11 = - facTU * 2./3. * real( conj(coupSUSYPtr->getLsqqG(iGen3,idIn1A))
1461                                        * coupSUSYPtr->getLsqqG(iGen4,idIn2A));
1462         GG22 = - facTU * 2./3. * real( conj(coupSUSYPtr->getRsqqG(iGen3,idIn1A))
1463                                        * coupSUSYPtr->getRsqqG(iGen4,idIn2A));
1464         // Sum over two contributing helicities
1465         sumInterference += sigmaGlu / sH / tGlu
1466           * ( GG11 + GG22 );
1467       }
1468
1469     }
1470
1471     // Skip the rest if only including QCD diagrams
1472     if (onlyQCD) return sumColT+sumColS+sumInterference;
1473
1474     // s-channel photon (strictly flavor-diagonal) and Z/gamma interference
1475     if (abs(id3Sav) == abs(id4Sav) && abs(id1) == abs(id2)) {
1476
1477       // gamma
1478       // Factor 2 since contributes to both ha != hb helicities
1479       sumColS += 2. * pow2(eQ) * pow2(eSq) * sigmaEW * facTU / pow2(sH);
1480
1481       // Z/gamma interference
1482       double CsqZ = real(coupSUSYPtr->LsusuZ[iGen3][iGen4]
1483                          + coupSUSYPtr->RsusuZ[iGen3][iGen4]);
1484       if (abs(id3Sav)%2 != 0) CsqZ = real(coupSUSYPtr->LsdsdZ[iGen3][iGen4]
1485                                           + coupSUSYPtr->RsdsdZ[iGen3][iGen4]);
1486       sumColS += eQ * eSq * sigmaEW * facTU / 2.0 / xW / (1.-xW)
1487         * sqrt(norm(propZW)) / sH * CsqZ
1488         * (coupSUSYPtr->LqqZ[idIn1A] + coupSUSYPtr->LqqZ[idIn2A]);
1489
1490       // Gluino/gamma interference (only for same-isospin)
1491       if (eQ == eSq) {
1492         double CsqG11 = real(conj(coupSUSYPtr->LsuuG[iGen3][iGen1])
1493                              *coupSUSYPtr->LsuuG[iGen4][iGen2]);
1494         double CsqG22 = real(conj(coupSUSYPtr->RsuuG[iGen3][iGen1])
1495                              *coupSUSYPtr->RsuuG[iGen4][iGen2]);
1496         if (id3Sav%2 != 0) {
1497           CsqG11 = real(conj(coupSUSYPtr->LsddG[iGen3][iGen1])
1498                         *coupSUSYPtr->LsddG[iGen4][iGen2]);
1499           CsqG22 = real(conj(coupSUSYPtr->RsddG[iGen3][iGen1])
1500                         *coupSUSYPtr->RsddG[iGen4][iGen2]);
1501         }
1502         sumColS += eQ * eSq * sigmaEWG * facTU
1503           * (CsqG11 + CsqG22) / sH / tGlu;
1504       }
1505     }
1506
1507     // s-channel Z (only for q flavor = qbar flavor)
1508     if (abs(id1) == abs(id2)) {
1509       double CsqZ = norm(coupSUSYPtr->LsusuZ[iGen3][iGen4]
1510                          + coupSUSYPtr->RsusuZ[iGen3][iGen4]);
1511       if (abs(id3Sav)%2 != 0) CsqZ = norm(coupSUSYPtr->LsdsdZ[iGen3][iGen4]
1512                                           + coupSUSYPtr->RsdsdZ[iGen3][iGen4]);
1513       sumColS += sigmaEW * facTU / 16.0 / pow2(xW) / pow2(1.0-xW)
1514         * norm(propZW) * CsqZ
1515         * ( pow2(coupSUSYPtr->LqqZ[idIn1A]) + pow2(coupSUSYPtr->RqqZ[idIn1A]) );
1516
1517       // Z/gluino interference (only for in-isospin = out-isospin)
1518       if (eQ == eSq) {
1519         double GZ11 = real(conj(coupSUSYPtr->getLsqqG(iGen3,idIn1A))
1520                            *coupSUSYPtr->getLsqqG(iGen4,idIn2A)
1521                            *(coupSUSYPtr->getLsqsqZ(id3Sav,id4Sav)
1522                              +coupSUSYPtr->getRsqsqZ(id3Sav,id4Sav)))
1523           *coupSUSYPtr->LqqZ[idIn1A];
1524         double GZ22 = real(conj(coupSUSYPtr->getRsqqG(iGen3,idIn1A))
1525                            *coupSUSYPtr->getRsqqG(iGen4,idIn2A)
1526                            *(coupSUSYPtr->getLsqsqZ(id3Sav,id4Sav)
1527                              +coupSUSYPtr->getRsqsqZ(id3Sav,id4Sav)))
1528           *coupSUSYPtr->RqqZ[idIn1A];
1529         sumColS += sigmaEWG * facTU / 4.0 / xW / (1.-xW)
1530           * ( GZ11 + GZ22 ) * sqrt(norm(propZW)) / tGlu;
1531       }
1532     }
1533
1534     // t-channel neutralinos
1535     // NOT YET IMPLEMENTED !
1536
1537   }
1538
1539   // Cross section
1540   double sigma = sumColS + sumColT + sumInterference;
1541
1542   // Return answer.
1543   return sigma;
1544
1545 }
1546
1547 //--------------------------------------------------------------------------
1548
1549 // Select identity, colour and anticolour.
1550
1551 void Sigma2qqbar2squarkantisquark::setIdColAcol() {
1552
1553   // Check if charge conjugate final state?
1554   isCC = false;
1555   if (isUD && ( (id1-1)%2 < 0 || (id2-1)%2 < 0 )) isCC = true;
1556
1557   //check if charge conjugate
1558   id3 = (isCC) ? -id3Sav : id3Sav;
1559   id4 = (isCC) ? -id4Sav : id4Sav;
1560
1561   // Set flavours.
1562   setId( id1, id2, id3, id4);
1563
1564   // Coded UD sigma is for udbar -> ~u~d'*. Swap t<->u for dbaru -> ~u~d'*.
1565   // Coded QQ sigma is for qqbar -> ~q~q*. Swap t<->u for qbarq -> ~q~q*.
1566   if (isUD) {
1567     swapTU = (abs(id1) % 2 != 0);
1568   } else {
1569     swapTU = (id1 < 0);
1570   }
1571
1572   // Select colour flow topology
1573   double R = rndmPtr->flat();
1574   double fracS = sumColS / (sumColS + sumColT) ;
1575   // S: color flow as in S-channel singlet
1576   if (R < fracS) {
1577     setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
1578     if (swapTU) setColAcol( 0, 1, 1, 0, 2, 0, 0, 2);
1579   }
1580   // T: color flow as in T-channel singlet
1581   else {
1582     setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
1583     if (swapTU) setColAcol( 0, 1, 2, 0, 2, 0, 0, 1);
1584   }
1585
1586   if (isCC) swapColAcol();
1587
1588 }
1589
1590 //==========================================================================
1591
1592 // Sigma2gg2squarkantisquark
1593 // Cross section for gg-initiated squark-antisquark production
1594
1595 //--------------------------------------------------------------------------
1596
1597 // Initialize process.
1598
1599 void Sigma2gg2squarkantisquark::initProc() {
1600
1601   //Typecast to the correct couplings
1602   coupSUSYPtr = (CoupSUSY*) couplingsPtr;
1603
1604   // Process Name
1605   nameSave = "g g -> "+particleDataPtr->name(abs(id3Sav))+" "
1606     +particleDataPtr->name(-abs(id4Sav));
1607
1608   // Squark pole mass
1609   m2Sq = pow2(particleDataPtr->m0(id3Sav));
1610
1611   // Secondary open width fraction.
1612   openFracPair = particleDataPtr->resOpenFrac(id3Sav, id4Sav);
1613
1614 }
1615
1616 //--------------------------------------------------------------------------
1617
1618 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
1619
1620 void Sigma2gg2squarkantisquark::sigmaKin() {
1621
1622   // Modified Mandelstam variables for massive kinematics with m3 = m4.
1623   // tHSq = tHat - m_squark^2; uHSq = uHat - m_squark^2.
1624   //  double s34Avg = 0.5 * (s3 + s4) - 0.25 * pow2(s3 - s4) / sH;
1625   double tHSq    = -0.5 * (sH - tH + uH);
1626   double uHSq    = -0.5 * (sH + tH - uH);
1627   // ! (NEED TO CHECK THAT THESE APPLIED CORRECTLY BELOW)   !
1628   // ! (PRELIMINARY CROSS-CHECKS WITH PYTHIA 6 COME OUT OK) !
1629
1630   // Helicity-independent prefactor
1631   double comFacHat = M_PI/sH2 * pow2(alpS) / 128.0
1632     * ( 24.0 * (1.0 - 2*tHSq*uHSq/sH2) - 8.0/3.0 );
1633
1634   // Helicity-dependent factors
1635   sigma = 0.0;
1636   for (int ha=-1;ha<=1;ha += 2) {
1637     for (int hb=-1;hb<=1;hb += 2) {
1638       // Divide by 4 for helicity average
1639       sigma += comFacHat / 4.0
1640         * ( (1.0-ha*hb)
1641             - 2.0 * sH*m2Sq/tHSq/uHSq
1642             * ( 1.0 - ha*hb - sH*m2Sq/tHSq/uHSq));
1643     }
1644   }
1645
1646 }
1647
1648 //--------------------------------------------------------------------------
1649
1650 // Select identity, colour and anticolour.
1651
1652 void Sigma2gg2squarkantisquark::setIdColAcol() {
1653
1654   // Set flavours.
1655   setId( id1, id2, id3Sav, id4Sav);
1656
1657   // Set color flow (random for now)
1658   double R = rndmPtr->flat();
1659   if (R < 0.5) setColAcol( 1, 2, 2, 3, 1, 0, 0, 3);
1660   else setColAcol( 1, 2, 3, 1, 3, 0, 0, 2);
1661
1662 }
1663
1664 //==========================================================================
1665
1666 // Sigma2qg2squarkgluino
1667 // Cross section for squark-gluino production
1668
1669 //--------------------------------------------------------------------------
1670
1671 // Initialize process.
1672
1673 void Sigma2qg2squarkgluino::initProc() {
1674
1675   //Typecast to the correct couplings
1676   coupSUSYPtr = (CoupSUSY*) couplingsPtr;
1677
1678   // Derive name
1679   nameSave = "q g -> "+particleDataPtr->name(abs(id3Sav))+" gluino + c.c.";
1680
1681   // Final-state mass squares
1682   m2Glu     = pow2(particleDataPtr->m0(1000021));
1683   m2Sq      = pow2(particleDataPtr->m0(id3Sav));
1684
1685   // Secondary open width fraction.
1686   openFracPair = particleDataPtr->resOpenFrac(id3Sav, 1000021);
1687
1688 }
1689
1690 //--------------------------------------------------------------------------
1691
1692 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
1693
1694 void Sigma2qg2squarkgluino::sigmaKin() {
1695
1696   // Common pre-factor
1697   comFacHat = (M_PI / sH2) * pow2(alpS) * 0.5 * openFracPair;
1698
1699   // Invariants (still with Pythia 6 sign convention)
1700   double tGlu = m2Glu-tH;
1701   double uGlu = m2Glu-uH;
1702   double tSq  = m2Sq-tH;
1703   double uSq  = m2Sq-uH;
1704
1705   // Color flow A: quark color annihilates with anticolor of g
1706   sigmaA = 0.5*4./9.* tGlu/sH + (tGlu*sH+2.*m2Glu*tSq)/pow2(tGlu) -
1707     ( (sH-m2Sq+m2Glu)*(-tSq)-sH*m2Glu )/sH/(-tGlu)
1708     + 0.5*1./2.*( tSq*(tH+2.*uH+m2Glu)-tGlu*(sH-2.*tSq)
1709                   + (-uGlu)*(tH+m2Glu+2.*m2Sq) )/2./tGlu/uSq;
1710   // Color flow B: quark and gluon colors iterchanged
1711   sigmaB =     4./9.*(-uGlu)*(uH+m2Sq)/pow2(uSq)
1712     + 1./18.* (sH*(uH+m2Glu) + 2.*(m2Sq-m2Glu)*uGlu)/sH/(-uSq)
1713     + 0.5*4./9.*tGlu/sH
1714     + 0.5*1./2.*(tSq*(tH+2.*uH+m2Glu)-tGlu*(sH-2.*tSq)
1715                  + (-uGlu)*(tH+m2Glu+2.*m2Sq))/2./tGlu/uSq;
1716
1717 }
1718
1719 double Sigma2qg2squarkgluino::sigmaHat() {
1720
1721   // Check whether right incoming flavor
1722   int idQA = (id1 == 21) ? abs(id2) : abs(id1);
1723   int idSqSM = id3Sav%1000000;
1724   if (idQA != idSqSM) return 0.0;
1725   // NOTE: ONLY WORKS FOR SLHA1 ENUMERATION !!!
1726   // (should replace this by squark mixing matrix squares)
1727
1728   return comFacHat * (sigmaA + sigmaB);
1729
1730 }
1731
1732 //--------------------------------------------------------------------------
1733
1734 // Select identity, colour and anticolour.
1735
1736 void Sigma2qg2squarkgluino::setIdColAcol() {
1737
1738   // Check if charge conjugate final state?
1739   int idQ = (id1 == 21) ? id2 : id1;
1740   id3 = (idQ > 0) ? id3Sav : -id3Sav;
1741   id4 = 1000021;
1742
1743   // Set flavors
1744   setId( id1, id2, id3, id4);
1745
1746   // Select color flow A or B (see above)
1747   double R = rndmPtr->flat()*(sigmaA+sigmaB);
1748   if (idQ == id1) {
1749     setColAcol(1,0,2,1,3,0,2,3);
1750     if (R > sigmaA) setColAcol(1,0,2,3,2,0,1,3);
1751   } else {
1752     setColAcol(2,1,1,0,3,0,2,3);
1753     if (R > sigmaB) setColAcol(2,3,1,0,2,0,1,3);
1754   }
1755   if (idQ < 0) swapColAcol();
1756
1757   // Use reflected kinematics if gq initial state
1758   if (id1 == 21) swapTU = true;
1759
1760 }
1761
1762 //==========================================================================
1763
1764 // Sigma2gg2gluinogluino
1765 // Cross section for gluino pair production from gg initial states
1766 // (validated against Pythia 6)
1767
1768 //--------------------------------------------------------------------------
1769
1770 // Initialize process.
1771
1772 void Sigma2gg2gluinogluino::initProc() {
1773
1774   //Typecast to the correct couplings
1775   coupSUSYPtr = (CoupSUSY*) couplingsPtr;
1776
1777   // Secondary open width fraction.
1778   openFracPair = particleDataPtr->resOpenFrac(1000021, 1000021);
1779
1780 }
1781
1782 //--------------------------------------------------------------------------
1783
1784 // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
1785
1786 void Sigma2gg2gluinogluino::sigmaKin() {
1787
1788   // Modified Mandelstam variables for massive kinematics with m3 = m4.
1789   // tHG = tHat - m_gluino^2; uHG = uHat - m_gluino^2.
1790   double s34Avg = 0.5 * (s3 + s4) - 0.25 * pow2(s3 - s4) / sH;
1791   double tHG    = -0.5 * (sH - tH + uH);
1792   double uHG    = -0.5 * (sH + tH - uH);
1793   double tHG2   = tHG * tHG;
1794   double uHG2   = uHG * uHG;
1795
1796   // Calculate kinematics dependence.
1797   sigTS  = (tHG * uHG - 2. * s34Avg * (tHG + 2. * s34Avg)) / tHG2
1798          + (tHG * uHG + s34Avg * (uHG - tHG)) / (sH * tHG);
1799   sigUS  = (tHG * uHG - 2. * s34Avg * (uHG + 2. * s34Avg)) / uHG2
1800          + (tHG * uHG + s34Avg * (tHG - uHG)) / (sH * uHG);
1801   sigTU  = 2. * tHG * uHG / sH2 + s34Avg * (sH - 4. * s34Avg)
1802          / (tHG * uHG);
1803   sigSum = sigTS + sigUS + sigTU;
1804
1805   // Answer contains factor 1/2 from identical gluinos.
1806   sigma  = (M_PI / sH2) * pow2(alpS) * (9./4.) * 0.5 * sigSum
1807          * openFracPair;
1808
1809 }
1810
1811 //--------------------------------------------------------------------------
1812
1813 // Select identity, colour and anticolour.
1814
1815 void Sigma2gg2gluinogluino::setIdColAcol() {
1816
1817   // Flavours are trivial.
1818   setId( id1, id2, 1000021, 1000021);
1819
1820   // Three colour flow topologies, each with two orientations.
1821   double sigRand = sigSum * rndmPtr->flat();
1822   if (sigRand < sigTS) setColAcol( 1, 2, 2, 3, 1, 4, 4, 3);
1823   else if (sigRand < sigTS + sigUS)
1824                        setColAcol( 1, 2, 3, 1, 3, 4, 4, 2);
1825   else                 setColAcol( 1, 2, 3, 4, 1, 4, 3, 2);
1826   if (rndmPtr->flat() > 0.5) swapColAcol();
1827
1828 }
1829
1830 //==========================================================================
1831
1832 // Sigma2qqbar2gluinogluino
1833 // Cross section for gluino pair production from qqbar initial states
1834 // (validated against Pythia 6)
1835
1836 //--------------------------------------------------------------------------
1837
1838 // Initialize process.
1839
1840 void Sigma2qqbar2gluinogluino::initProc() {
1841
1842   //Typecast to the correct couplings
1843   coupSUSYPtr = (CoupSUSY*) couplingsPtr;
1844
1845   // Secondary open width fraction.
1846   openFracPair = particleDataPtr->resOpenFrac(1000021, 1000021);
1847
1848 }
1849
1850 //--------------------------------------------------------------------------
1851
1852 // Begin evaluate d(sigmaHat)/d(tHat); flavour-independent part.
1853
1854 void Sigma2qqbar2gluinogluino::sigmaKin() {
1855
1856   // Modified Mandelstam variables for massive kinematics with m3 = m4.
1857   // tHG = tHat - m_gluino^2; uHG = uHat - m_gluino^2.
1858   // (Note: tHG and uHG defined with opposite sign wrt Pythia 6)
1859   s34Avg = 0.5 * (s3 + s4) - 0.25 * pow2(s3 - s4) / sH;
1860   tHG    = -0.5 * (sH - tH + uH);
1861   uHG    = -0.5 * (sH + tH - uH);
1862   tHG2   = tHG * tHG;
1863   uHG2   = uHG * uHG;
1864
1865   // s-channel contribution.
1866   sigS   = (tHG2 + uHG2 + 2. * s34Avg * sH) / sH2;
1867
1868 }
1869
1870 //--------------------------------------------------------------------------
1871
1872
1873 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1874
1875 double Sigma2qqbar2gluinogluino::sigmaHat() {
1876
1877   // Squarks (L/R or 1/2) can contribute in t or u channel.
1878   // Assume identical CKM matrices in quark and squark sector.
1879   // (Note: tHQL,R and uHQL,R defined with opposite sign wrt Pythia 6.
1880   //  This affects the sign of the interference term below)
1881   double sQL    = pow2( particleDataPtr->m0(1000000 + abs(id1)) );
1882   double tHQL   = tHG + s34Avg - sQL;
1883   double uHQL   = uHG + s34Avg - sQL;
1884   double sQR    = pow2( particleDataPtr->m0(2000000 + abs(id1)) );
1885   double tHQR   = tHG + s34Avg - sQR;
1886   double uHQR   = uHG + s34Avg - sQR;
1887
1888   // Calculate kinematics dependence.
1889   double sigQL  = (4./9.) * (tHG2 / pow2(tHQL) + uHG2 / pow2(uHQL))
1890                 + (1./9.) * s34Avg * sH / (tHQL * uHQL)
1891                 + (tHG2 + sH * s34Avg) /(sH * tHQL)
1892                 + (uHG2 + sH * s34Avg) /(sH * uHQL);
1893   double sigQR  = (4./9.) * (tHG2 / pow2(tHQR) + uHG2 / pow2(uHQR))
1894                 + (1./9.) * s34Avg * sH / (tHQR * uHQR)
1895                 + (tHG2 + sH * s34Avg) /(sH * tHQR)
1896                 + (uHG2 + sH * s34Avg) /(sH * uHQR);
1897   double sigSum = sigS + 0.5 * (sigQL + sigQR);
1898
1899   // Answer contains factor 1/2 from identical gluinos.
1900   double sigma  = (M_PI / sH2) * pow2(alpS) * (8./3.) * 0.5 * sigSum
1901                 * openFracPair;
1902   return sigma;
1903
1904 }
1905
1906 //--------------------------------------------------------------------------
1907
1908 // Select identity, colour and anticolour.
1909
1910 void Sigma2qqbar2gluinogluino::setIdColAcol() {
1911
1912   // Flavours are trivial.
1913   setId( id1, id2, 1000021, 1000021);
1914
1915   // Two colour flow topologies. Swap if first is antiquark.
1916   if (rndmPtr->flat() < 0.5) setColAcol( 1, 0, 0, 2, 1, 3, 3, 2);
1917   else                    setColAcol( 1, 0, 0, 2, 3, 2, 1, 3);
1918   if (id1 < 0) swapColAcol();
1919
1920 }
1921
1922 //==========================================================================
1923
1924 // Sigma1qq2antisquark
1925 // R-parity violating squark production
1926
1927 //--------------------------------------------------------------------------
1928
1929 // Initialise process
1930
1931 void Sigma1qq2antisquark::initProc(){
1932
1933   //Typecast to the correct couplings
1934   coupSUSYPtr = (CoupSUSY*) couplingsPtr;
1935
1936   //Construct name of the process from lambda'' couplings
1937
1938   nameSave = "q q' -> " + coupSUSYPtr->getName(idRes)+"* + c.c";
1939   codeSave = 2000 + 10*abs(idRes)/1000000 + abs(idRes)%10;
1940 }
1941
1942 //--------------------------------------------------------------------------
1943
1944 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
1945
1946 void Sigma1qq2antisquark::sigmaKin() {
1947
1948   // Check if at least one RPV coupling non-zero
1949   if(!coupSUSYPtr->isUDD) {
1950     sigBW = 0.0;
1951     return;
1952   }
1953
1954   mRes = particleDataPtr->m0(abs(idRes));
1955   GammaRes = particleDataPtr->mWidth(abs(idRes));
1956   m2Res = pow2(mRes);
1957
1958   sigBW        = sH * GammaRes/ ( pow2(sH - m2Res) + pow2(mRes * GammaRes) );
1959   sigBW       *= 2.0/3.0/mRes;
1960
1961   // Width out only includes open channels.
1962   widthOut     = GammaRes * particleDataPtr->resOpenFrac(id3);
1963 }
1964
1965 //--------------------------------------------------------------------------
1966
1967 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1968
1969 double Sigma1qq2antisquark::sigmaHat() {
1970
1971   // Only allow (anti)quark-(anti)quark incoming states
1972   if (id1*id2 <= 0) return 0.0;
1973
1974   // Generation indices
1975   int iA = (abs(id1)+1)/2;
1976   int iB = (abs(id2)+1)/2;
1977
1978   //Covert from pdg-code to ~u_i/~d_i basis
1979   bool idown = (abs(idRes)%2 == 1) ? true : false;
1980   int iC = (abs(idRes)/1000000 == 2) ? (abs(idRes)%10+1)/2 + 3: (abs(idRes)%10+1)/2;
1981
1982   // UDD structure
1983   if (abs(id1)%2 == 0 && abs(id2)%2 == 0) return 0.0;
1984   if (abs(id1)%2 == 1 && abs(id2)%2 == 1 && idown) return 0.0;
1985   if ((abs(id1) + abs(id2))%2 == 1 && !idown) return 0.0;
1986
1987   double sigma = 0.0;
1988
1989   if(!idown){
1990    // d_i d_j --> ~u*_k
1991     for(int isq = 1; isq <=3; isq++){
1992       // Loop over R-type squark contributions
1993       sigma += pow2(coupSUSYPtr->rvUDD[isq][iA][iB])
1994         * norm(coupSUSYPtr->Rusq[iC][isq+3]);
1995     }
1996   }else{
1997     // u_i d_j --> d*_k
1998     // Pick the correct coupling for d-u in-state
1999     if(abs(id1)%2==1){
2000       iA = (abs(id2)+1)/2;
2001       iB = (abs(id1)+1)/2;
2002     }
2003     for(int isq = 1; isq <= 3; isq++){
2004       // Loop over R-type squark contributions
2005       sigma += pow2(coupSUSYPtr->rvUDD[iA][iB][isq])
2006         * norm(coupSUSYPtr->Rdsq[iC][isq+3]);
2007     }
2008   }
2009
2010   sigma *= sigBW;
2011   return sigma;
2012
2013 }
2014
2015 //--------------------------------------------------------------------------
2016
2017 // Select identity, colour and anticolour.
2018
2019 void Sigma1qq2antisquark::setIdColAcol() {
2020
2021   // Set flavours.
2022   if(id1 < 0 && id2 < 0 ) setId( id1, id2, idRes);
2023   else setId( id1, id2, -idRes);
2024
2025   // Colour flow topologies. Swap when antiquarks.
2026   if (abs(id1) < 9) setColAcol( 1, 0, 2, 0, 0, 3);
2027   else              setColAcol( 0, 0, 0, 0, 0, 0);
2028   if (id1 < 0) swapColAcol();
2029
2030 }
2031
2032
2033 //==========================================================================
2034
2035 } // end namespace Pythia8
2036