PYTHIA8/pythia8170/src/SigmaSUSY.cxx

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