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