]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PYTHIA8/pythia8170/src/SigmaSUSY.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Stupid bug fix in new superlight mode (from Zurich airport)
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8170 / src / SigmaSUSY.cxx
CommitLineData
63ba5337 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
12namespace Pythia8 {
13
14//==========================================================================
15
16// Sigma2qqbar2chi0chi0
17// Cross section for gaugino pair production: neutralino pair
18
19//--------------------------------------------------------------------------
20
21// Initialize process.
22
23void 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
41void 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
62double 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
179void 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
200void 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
222double 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
329void 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
349double 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
483void 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
507void 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
526double 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
594void 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
615void 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
639double 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
705void 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
730void 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
774void 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
803double 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
1255void 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
1291void 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
1334void 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
1360double 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
1567void 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
1615void 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
1636void 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
1668void 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
1689void 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
1710void 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
1735double 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
1752void 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
1788void 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
1802void 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
1831void 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
1856void 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
1870void 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
1891double 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
1926void 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
1947void 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
1962void 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
1985double 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
2036void 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