Line data Source code
1 : // SigmaSUSY.cc is a part of the PYTHIA event generator.
2 : // Copyright (C) 2015 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 "Pythia8/SigmaSUSY.h"
11 :
12 : namespace Pythia8 {
13 :
14 : //==========================================================================
15 :
16 : // Sigma2SUSY
17 : // An intermediate class for SUSY 2 -> 2 with nontrivial decay angles.
18 :
19 : //--------------------------------------------------------------------------
20 :
21 : double Sigma2SUSY::weightDecay( Event& process, int iResBeg, int iResEnd) {
22 :
23 : // Do nothing if decays present already at input.
24 :
25 : // Identity of mother of decaying resonance(s).
26 0 : int idMother = process[process[iResBeg].mother1()].idAbs();
27 :
28 : // For Higgs decay hand over to standard routine.
29 0 : if (idMother == 25 || idMother == 35 || idMother == 36)
30 0 : return weightHiggsDecay( process, iResBeg, iResEnd);
31 :
32 : // For top decay hand over to standard routine.
33 0 : if (idMother == 6)
34 0 : return weightTopDecay( process, iResBeg, iResEnd);
35 :
36 : // For Neutralino(i) decay hand over to standard routine.
37 0 : if ( settingsPtr->flag("SUSYResonance:3BodyMatrixElement")
38 0 : && (idMother == 1000023 || idMother == 1000025 || idMother == 1000035) ) {
39 :
40 : // Nj -> Ni f fbar
41 0 : if (iResEnd - iResBeg != 2) return(1.0);
42 : int iW1 = iResBeg;
43 0 : int iF = iResBeg + 1;
44 0 : int iFbar= iResBeg + 2;
45 0 : int iT = process[iW1].mother1();
46 0 : if( iT <= 0 ) return(1.0);
47 0 : int idDau= process[iW1].idAbs();
48 :
49 : // Neutralino decays to charginos not yet implemented
50 0 : if (idDau == 1000024 || idDau == 1000037 ) return(1.0);
51 0 : if (idDau != 1000022 && idDau != 1000023 && idDau != 1000025
52 0 : && idDau != 1000035 ) {
53 0 : return(1.0);
54 : } else {
55 0 : if( process[iF].idAbs() != process[iFbar].idAbs() ) return(1.0);
56 : int idmo = -1; int iddau = -1;
57 0 : switch (idMother) {
58 0 : case 1000023: idmo = 2; break;
59 0 : case 1000025: idmo = 3; break;
60 0 : case 1000035: idmo = 4; break;
61 : }
62 0 : switch (idDau) {
63 0 : case 1000022: iddau = 1; break;
64 0 : case 1000023: iddau = 2; break;
65 0 : case 1000025: iddau = 3; break;
66 : }
67 0 : if( idmo<0 || iddau<0 ) return(1.0);
68 :
69 0 : Sigma2qqbar2chi0chi0 localDecay(idmo,iddau,0);
70 0 : localDecay.init(infoPtr, settingsPtr, particleDataPtr,NULL,NULL,
71 0 : NULL,couplingsPtr);
72 0 : localDecay.initProc();
73 0 : localDecay.alpEM = 1;
74 0 : localDecay.id1 = process[iF].id();
75 0 : localDecay.id2 = process[iFbar].id();
76 0 : double xm3 = process[iT].m();
77 0 : double xm4 = process[iW1].m();
78 0 : localDecay.m3 = xm3;
79 0 : localDecay.m4 = xm4;
80 0 : localDecay.s3 = xm3*xm3;
81 0 : localDecay.s4 = xm4*xm4;
82 0 : localDecay.sH = (process[iF].p()+process[iFbar].p()).m2Calc();
83 0 : localDecay.sH2 = pow2(localDecay.sH);
84 0 : localDecay.tH = (process[iF].p()-process[iT].p()).m2Calc();
85 0 : localDecay.uH = localDecay.s3+localDecay.s4-localDecay.tH-localDecay.sH;
86 0 : localDecay.sigmaKin();
87 0 : double wt = -localDecay.sigmaHat();
88 : // Estimate maximum weight by sampling kinematic extremes
89 : // Case I: neutralino(i) at rest
90 0 : localDecay.sH = pow2(xm4-xm3);
91 0 : localDecay.tH = 0.5*(localDecay.s3+localDecay.s4-localDecay.sH);
92 0 : localDecay.uH = localDecay.tH;
93 0 : localDecay.sigmaKin();
94 0 : double wtmax = -localDecay.sigmaHat();
95 : // Case II: fermion at rest
96 0 : localDecay.sH = 0;
97 0 : localDecay.tH = localDecay.s3;
98 0 : localDecay.uH = localDecay.s3+localDecay.s4-localDecay.tH-localDecay.sH;
99 0 : localDecay.sigmaKin();
100 0 : wtmax += -localDecay.sigmaHat();
101 : // Case III: antifermion at rest
102 0 : localDecay.uH = localDecay.s3;
103 0 : localDecay.tH = localDecay.s3+localDecay.s4-localDecay.tH-localDecay.sH;
104 0 : localDecay.sigmaKin();
105 0 : wtmax += -localDecay.sigmaHat();
106 0 : return(wt/wtmax);
107 0 : }
108 : }
109 :
110 : // Else done.
111 0 : return 1.;
112 :
113 0 : }
114 :
115 : //==========================================================================
116 :
117 : // Sigma2qqbar2chi0chi0
118 : // Cross section for gaugino pair production: neutralino pair
119 :
120 : //--------------------------------------------------------------------------
121 :
122 : // Initialize process.
123 :
124 : void Sigma2qqbar2chi0chi0::initProc() {
125 :
126 : //Typecast to the correct couplings
127 0 : coupSUSYPtr = (CoupSUSY*) couplingsPtr;
128 :
129 : // Construct name of process.
130 0 : nameSave = "q qbar' -> " + particleDataPtr->name(id3) + " "
131 0 : + particleDataPtr->name(id4);
132 :
133 : // Secondary open width fraction.
134 0 : openFracPair = particleDataPtr->resOpenFrac(id3, id4);
135 :
136 0 : }
137 :
138 : //--------------------------------------------------------------------------
139 :
140 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
141 :
142 : void Sigma2qqbar2chi0chi0::sigmaKin() {
143 :
144 : // Common flavour-independent factor.
145 0 : sigma0 = M_PI /3.0/ sH2 / pow2(coupSUSYPtr->sin2W) * pow2(alpEM)
146 0 : * openFracPair;
147 :
148 : // Auxiliary factors for use below
149 0 : ui = uH - s3;
150 0 : uj = uH - s4;
151 0 : ti = tH - s3;
152 0 : tj = tH - s4;
153 0 : double sV= sH - pow2(coupSUSYPtr->mZpole);
154 0 : double d = pow2(sV) + pow2(coupSUSYPtr->mZpole * coupSUSYPtr->wZpole);
155 0 : propZ = complex( sV / d, coupSUSYPtr->mZpole * coupSUSYPtr->wZpole / d);
156 :
157 0 : }
158 :
159 : //--------------------------------------------------------------------------
160 :
161 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
162 :
163 : double Sigma2qqbar2chi0chi0::sigmaHat() {
164 :
165 : // Only allow quark-antiquark incoming states
166 0 : if (id1*id2 >= 0) {
167 0 : return 0.0;
168 : }
169 :
170 : // Only allow incoming states with sum(charge) = 0
171 0 : if ((id1+id2) % 2 != 0) {
172 0 : return 0.0;
173 : }
174 :
175 0 : if(id1<0) swapTU = true;
176 :
177 : // Shorthands
178 0 : int idAbs1 = abs(id1);
179 0 : int idAbs2 = abs(id2);
180 :
181 : // Flavour-dependent kinematics-dependent couplings.
182 0 : complex QuLL(0.0),QtLL(0.0),QuRR(0.0),QtRR(0.0);
183 0 : complex QuLR(0.0),QtLR(0.0),QuRL(0.0),QtRL(0.0);
184 :
185 : double *LqqZloc;
186 : double *RqqZloc;
187 :
188 : int iAdd=0;
189 :
190 0 : if( idAbs1 > 10 && idAbs1 < 17 ) {
191 0 : LqqZloc = coupSUSYPtr->LllZ;
192 0 : RqqZloc = coupSUSYPtr->RllZ;
193 : iAdd+=10;
194 0 : } else {
195 0 : LqqZloc = coupSUSYPtr->LqqZ;
196 0 : RqqZloc = coupSUSYPtr->RqqZ;
197 : }
198 :
199 : // s-channel Z couplings
200 0 : if (idAbs1 == idAbs2) {
201 0 : QuLL = LqqZloc[idAbs1-iAdd] * coupSUSYPtr->OLpp[id3chi][id4chi]
202 0 : * propZ / 2.0;
203 0 : QtLL = LqqZloc[idAbs1-iAdd] * coupSUSYPtr->ORpp[id3chi][id4chi]
204 0 : * propZ / 2.0;
205 0 : QuRR = RqqZloc[idAbs1-iAdd] * coupSUSYPtr->ORpp[id3chi][id4chi]
206 0 : * propZ / 2.0;
207 0 : QtRR = RqqZloc[idAbs1-iAdd] * coupSUSYPtr->OLpp[id3chi][id4chi]
208 0 : * propZ / 2.0;
209 0 : }
210 :
211 : // Flavour indices
212 0 : int ifl1 = (idAbs1+1-iAdd) / 2;
213 0 : int ifl2 = (idAbs2+1-iAdd) / 2;
214 :
215 : complex (*LsddXloc)[4][6];
216 : complex (*RsddXloc)[4][6];
217 : complex (*LsuuXloc)[4][6];
218 : complex (*RsuuXloc)[4][6];
219 0 : if( idAbs1 > 10 && idAbs1 < 17 ) {
220 0 : LsddXloc = coupSUSYPtr->LsllX;
221 0 : RsddXloc = coupSUSYPtr->RsllX;
222 0 : LsuuXloc = coupSUSYPtr->LsvvX;
223 0 : RsuuXloc = coupSUSYPtr->RsvvX;
224 0 : } else {
225 0 : LsddXloc = coupSUSYPtr->LsddX;
226 0 : RsddXloc = coupSUSYPtr->RsddX;
227 0 : LsuuXloc = coupSUSYPtr->LsuuX;
228 0 : RsuuXloc = coupSUSYPtr->RsuuX;
229 : }
230 :
231 : // Add t-channel squark flavour sums to QmXY couplings
232 0 : for (int ksq=1; ksq<=6; ksq++) {
233 :
234 : // squark id and squark-subtracted u and t
235 :
236 : int idsq;
237 0 : idsq=((ksq+2)/3)*1000000 + 2*((ksq-1) % 3) + (idAbs1+1) % 2 + 1;
238 0 : idsq+=iAdd;
239 :
240 0 : double msq2 = pow(particleDataPtr->m0(idsq),2);
241 0 : double usq = uH - msq2;
242 0 : double tsq = tH - msq2;
243 :
244 0 : complex Lsqq1X3;
245 0 : complex Lsqq1X4;
246 0 : complex Lsqq2X3;
247 0 : complex Lsqq2X4;
248 0 : complex Rsqq1X3;
249 0 : complex Rsqq1X4;
250 0 : complex Rsqq2X3;
251 0 : complex Rsqq2X4;
252 :
253 : // Couplings
254 0 : Lsqq1X3 = LsuuXloc[ksq][ifl1][id3chi];
255 0 : Lsqq1X4 = LsuuXloc[ksq][ifl1][id4chi];
256 0 : Lsqq2X3 = LsuuXloc[ksq][ifl2][id3chi];
257 0 : Lsqq2X4 = LsuuXloc[ksq][ifl2][id4chi];
258 0 : Rsqq1X3 = RsuuXloc[ksq][ifl1][id3chi];
259 0 : Rsqq1X4 = RsuuXloc[ksq][ifl1][id4chi];
260 0 : Rsqq2X3 = RsuuXloc[ksq][ifl2][id3chi];
261 0 : Rsqq2X4 = RsuuXloc[ksq][ifl2][id4chi];
262 0 : if (idAbs1 % 2 != 0) {
263 0 : Lsqq1X3 = LsddXloc[ksq][ifl1][id3chi];
264 0 : Lsqq1X4 = LsddXloc[ksq][ifl1][id4chi];
265 0 : Lsqq2X3 = LsddXloc[ksq][ifl2][id3chi];
266 0 : Lsqq2X4 = LsddXloc[ksq][ifl2][id4chi];
267 0 : Rsqq1X3 = RsddXloc[ksq][ifl1][id3chi];
268 0 : Rsqq1X4 = RsddXloc[ksq][ifl1][id4chi];
269 0 : Rsqq2X3 = RsddXloc[ksq][ifl2][id3chi];
270 0 : Rsqq2X4 = RsddXloc[ksq][ifl2][id4chi];
271 0 : }
272 :
273 : // QuXY
274 0 : QuLL += conj(Lsqq1X4)*Lsqq2X3/usq;
275 0 : QuRR += conj(Rsqq1X4)*Rsqq2X3/usq;
276 0 : QuLR += conj(Lsqq1X4)*Rsqq2X3/usq;
277 0 : QuRL += conj(Rsqq1X4)*Lsqq2X3/usq;
278 :
279 : // QtXY
280 0 : QtLL -= conj(Lsqq1X3)*Lsqq2X4/tsq;
281 0 : QtRR -= conj(Rsqq1X3)*Rsqq2X4/tsq;
282 0 : QtLR += conj(Lsqq1X3)*Rsqq2X4/tsq;
283 0 : QtRL += conj(Rsqq1X3)*Lsqq2X4/tsq;
284 :
285 0 : }
286 :
287 : // Overall factor multiplying each coupling; multiplied at the end as fac^2
288 0 : double fac = (1.0-coupSUSYPtr->sin2W);
289 0 : if(abs(id3)==abs(id4)) fac *= sqrt(2.); // for identical final particles
290 :
291 : // Compute matrix element weight
292 : double weight = 0;
293 0 : double facLR = uH*tH - s3*s4;
294 0 : double facMS = m3*m4*sH;
295 :
296 : // Average over separate helicity contributions
297 : // LL (ha = -1, hb = +1) (divided by 4 for average)
298 0 : weight += norm(QuLL) * ui * uj + norm(QtLL) * ti * tj
299 0 : + 2 * real(conj(QuLL) * QtLL) * facMS;
300 : // RR (ha = 1, hb = -1) (divided by 4 for average)
301 0 : weight += norm(QtRR) * ti * tj + norm(QuRR) * ui * uj
302 0 : + 2 * real(conj(QuRR) * QtRR) * facMS;
303 : // RL (ha = 1, hb = 1) (divided by 4 for average)
304 0 : weight += norm(QuRL) * ui * uj + norm(QtRL) * ti * tj
305 0 : + real(conj(QuRL) * QtRL) * facLR;
306 : // LR (ha = -1, hb = -1) (divided by 4 for average)
307 0 : weight += norm(QuLR) * ui * uj + norm(QtLR) * ti * tj
308 0 : + real(conj(QuLR) * QtLR) * facLR;
309 :
310 0 : double colorFactor = ( idAbs1 > 10 && idAbs1 < 17 ) ? 3.0 : 1.0;
311 :
312 : // Cross section, including colour factor.
313 0 : double sigma = sigma0 * weight / pow2(fac) * colorFactor;
314 :
315 : // Answer.
316 : return sigma;
317 :
318 0 : }
319 :
320 : //--------------------------------------------------------------------------
321 :
322 : // Select identity, colour and anticolour.
323 :
324 : void Sigma2qqbar2chi0chi0::setIdColAcol() {
325 :
326 : // Set flavours.
327 0 : setId( id1, id2, id3, id4);
328 :
329 : // Colour flow topologies. Swap when antiquarks.
330 0 : if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
331 0 : else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
332 0 : if (id1 < 0) swapColAcol();
333 :
334 0 : }
335 :
336 : //==========================================================================
337 :
338 : // Sigma2qqbar2charchi0
339 : // Cross section for gaugino pair production: neutralino-chargino
340 :
341 : //--------------------------------------------------------------------------
342 :
343 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
344 :
345 : void Sigma2qqbar2charchi0::sigmaKin() {
346 :
347 : // Common flavour-independent factor.
348 :
349 0 : sigma0 = M_PI / sH2 / 3.0 / pow2(coupSUSYPtr->sin2W) * pow2(alpEM) ;
350 0 : sigma0 /= 2.0 * (1 - coupSUSYPtr->sin2W) ;
351 :
352 : // Auxiliary factors for use below
353 0 : ui = uH - s3;
354 0 : uj = uH - s4;
355 0 : ti = tH - s3;
356 0 : tj = tH - s4;
357 0 : double sW = sH - pow2(coupSUSYPtr->mWpole);
358 0 : double d = pow2(sW) + pow2(coupSUSYPtr->mWpole * coupSUSYPtr->wWpole);
359 0 : propW = complex( sW / d, coupSUSYPtr->mWpole * coupSUSYPtr->wWpole / d);
360 :
361 0 : }
362 :
363 : //--------------------------------------------------------------------------
364 :
365 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
366 :
367 : double Sigma2qqbar2charchi0::sigmaHat() {
368 :
369 : // Only allow particle-antiparticle incoming states
370 0 : if (id1*id2 >= 0) {
371 0 : return 0.0;
372 : }
373 :
374 : // Only allow incoming states with sum(charge) = final state
375 0 : if (abs(id1) % 2 == abs(id2) % 2) return 0.0;
376 0 : int isPos = (id3chi > 0 ? 1 : 0);
377 0 : if (id1 < 0 && id1 > -19 && abs(id1) % 2 == 1-isPos ) return 0.0;
378 0 : else if (id1 > 0 && id1 < 19 && abs(id1) % 2 == isPos ) return 0.0;
379 :
380 : // Flavour-dependent kinematics-dependent couplings.
381 0 : int idAbs1 = abs(id1);
382 0 : int iChar = abs(id3chi);
383 0 : int iNeut = abs(id4chi);
384 :
385 0 : complex QuLL(0.0),QtLL(0.0),QuRR(0.0),QtRR(0.0);
386 0 : complex QuLR(0.0),QtLR(0.0),QuRL(0.0),QtRL(0.0);
387 :
388 : // Calculate everything from udbar -> ~chi+ ~chi0 template process
389 : int iAdd=0;
390 : complex (*LudWloc)[4];
391 : complex (*LsddXloc)[4][6];
392 : complex (*RsddXloc)[4][6];
393 : complex (*LsuuXloc)[4][6];
394 : complex (*RsuuXloc)[4][6];
395 : complex (*LsduXloc)[4][3];
396 : complex (*RsduXloc)[4][3];
397 : complex (*LsudXloc)[4][3];
398 : complex (*RsudXloc)[4][3];
399 0 : if( idAbs1 > 10 && idAbs1 < 17 ) {
400 : iAdd+=10;
401 0 : LudWloc = coupSUSYPtr->LlvW;
402 0 : LsddXloc = coupSUSYPtr->LsllX;
403 0 : RsddXloc = coupSUSYPtr->RsllX;
404 0 : LsuuXloc = coupSUSYPtr->LsvvX;
405 0 : RsuuXloc = coupSUSYPtr->RsvvX;
406 0 : LsduXloc = coupSUSYPtr->LslvX;
407 0 : RsduXloc = coupSUSYPtr->RslvX;
408 0 : LsudXloc = coupSUSYPtr->LsvlX;
409 0 : RsudXloc = coupSUSYPtr->RsvlX;
410 0 : } else {
411 0 : LudWloc = coupSUSYPtr->LudW;
412 0 : LsddXloc = coupSUSYPtr->LsddX;
413 0 : RsddXloc = coupSUSYPtr->RsddX;
414 0 : LsuuXloc = coupSUSYPtr->LsuuX;
415 0 : RsuuXloc = coupSUSYPtr->RsuuX;
416 0 : LsduXloc = coupSUSYPtr->LsduX;
417 0 : RsduXloc = coupSUSYPtr->RsduX;
418 0 : LsudXloc = coupSUSYPtr->LsudX;
419 0 : RsudXloc = coupSUSYPtr->RsudX;
420 : }
421 :
422 : // u dbar , ubar d : do nothing
423 : // dbar u , d ubar : swap 1<->2 and t<->u
424 0 : int iGu = (abs(id1)-iAdd)/2;
425 0 : int iGd = (abs(id2)+1-iAdd)/2;
426 0 : if (idAbs1 % 2 != 0) {
427 0 : swapTU = true;
428 0 : iGu = (abs(id2)-iAdd)/2;
429 0 : iGd = (abs(id1)+1-iAdd)/2;
430 0 : }
431 :
432 : // s-channel W contribution
433 0 : QuLL = conj(LudWloc[iGu][iGd])
434 0 : * conj(coupSUSYPtr->OL[iNeut][iChar])
435 0 : * propW / sqrt(2.0);
436 0 : QtLL = conj(LudWloc[iGu][iGd])
437 0 : * conj(coupSUSYPtr->OR[iNeut][iChar])
438 0 : * propW / sqrt(2.0);
439 :
440 : // Add t-channel squark flavour sums to QmXY couplings
441 0 : for (int jsq=1; jsq<=6; jsq++) {
442 :
443 0 : int idsu=((jsq+2)/3)*1000000 + 2*((jsq-1) % 3) + 2 +iAdd;
444 0 : int idsd=((jsq+2)/3)*1000000 + 2*((jsq-1) % 3) + 1 +iAdd;
445 0 : double msd2 = pow(particleDataPtr->m0(idsd),2);
446 0 : double msu2 = pow(particleDataPtr->m0(idsu),2);
447 0 : double tsq = tH - msd2;
448 0 : double usq = uH - msu2;
449 :
450 0 : QuLL += conj(LsuuXloc[jsq][iGu][iNeut])
451 0 : * conj(LsudXloc[jsq][iGd][iChar])/usq;
452 0 : QuLR += conj(LsuuXloc[jsq][iGu][iNeut])
453 0 : * conj(RsudXloc[jsq][iGd][iChar])/usq;
454 0 : QuRR += conj(RsuuXloc[jsq][iGu][iNeut])
455 0 : * conj(RsudXloc[jsq][iGd][iChar])/usq;
456 0 : QuRL += conj(RsuuXloc[jsq][iGu][iNeut])
457 0 : * conj(LsudXloc[jsq][iGd][iChar])/usq;
458 :
459 0 : QtLL -= conj(LsduXloc[jsq][iGu][iChar])
460 0 : * LsddXloc[jsq][iGd][iNeut]/tsq;
461 0 : QtRR -= conj(RsduXloc[jsq][iGu][iChar])
462 0 : * RsddXloc[jsq][iGd][iNeut]/tsq;
463 0 : QtLR += conj(LsduXloc[jsq][iGu][iChar])
464 0 : * RsddXloc[jsq][iGd][iNeut]/tsq;
465 0 : QtRL += conj(RsduXloc[jsq][iGu][iChar])
466 0 : * LsddXloc[jsq][iGd][iNeut]/tsq;
467 0 : }
468 :
469 : // Compute matrix element weight
470 : double weight = 0;
471 :
472 : // Average over separate helicity contributions
473 : // (if swapped, swap ha, hb if computing polarized cross sections)
474 : // LL (ha = -1, hb = +1) (divided by 4 for average)
475 0 : weight += norm(QuLL) * ui * uj + norm(QtLL) * ti * tj
476 0 : + 2 * real(conj(QuLL) * QtLL) * m3 * m4 * sH;
477 : // RR (ha = 1, hb = -1) (divided by 4 for average)
478 0 : weight += norm(QtRR) * ti * tj + norm(QuRR) * ui * uj
479 0 : + 2 * real(conj(QuRR) * QtRR) * m3 * m4 * sH;
480 : // RL (ha = 1, hb = 1) (divided by 4 for average)
481 0 : weight += norm(QuRL) * ui * uj + norm(QtRL) * ti * tj
482 0 : + real(conj(QuRL) * QtRL) * (uH * tH - s3 * s4);
483 : // LR (ha = -1, hb = -1) (divided by 4 for average)
484 0 : weight += norm(QuLR) * ui * uj + norm(QtLR) * ti * tj
485 0 : + real(conj(QuLR) * QtLR) * (uH * tH - s3 * s4);
486 :
487 0 : double colorFactor = ( idAbs1 > 10 && idAbs1 < 17 ) ? 3.0 : 1.0;
488 :
489 : // Cross section, including colour factor.
490 0 : double sigma = sigma0 * weight * colorFactor;
491 :
492 : // Answer.
493 : return sigma;
494 :
495 0 : }
496 :
497 : //==========================================================================
498 :
499 : // Sigma2qqbar2charchar
500 : // Cross section for gaugino pair production: chargino-chargino
501 :
502 : //--------------------------------------------------------------------------
503 :
504 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
505 :
506 : void Sigma2qqbar2charchar::sigmaKin() {
507 :
508 : // Common flavour-independent factor.
509 0 : sigma0 = M_PI / 3.0 / sH2 / pow2(coupSUSYPtr->sin2W) * pow2(alpEM) ;
510 :
511 : // Auxiliary factors for use below
512 0 : ui = uH - s3;
513 0 : uj = uH - s4;
514 0 : ti = tH - s3;
515 0 : tj = tH - s4;
516 0 : double sV= sH - pow2(coupSUSYPtr->mZpole);
517 0 : double d = pow2(sV) + pow2(coupSUSYPtr->mZpole * coupSUSYPtr->wZpole);
518 0 : propZ = complex( sV / d, coupSUSYPtr->mZpole * coupSUSYPtr->wZpole / d);
519 :
520 0 : }
521 :
522 : //--------------------------------------------------------------------------
523 :
524 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
525 :
526 : double Sigma2qqbar2charchar::sigmaHat() {
527 :
528 : // Only allow quark-antiquark incoming states
529 0 : if (id1*id2 >= 0) return 0.0;
530 :
531 : // Only allow incoming states with sum(charge) = 0
532 0 : if ((id1+id2) % 2 != 0) return 0.0;
533 :
534 : //if (id1 > 0 || id1==-1 || id1==-3 || id1==-5) return 0.0;
535 : //if (id1 < 0 || id1==1 || id1==3 || id1==5) return 0.0;
536 :
537 0 : swapTU = (id1 < 0 ? true : false);
538 :
539 : // Flavour-dependent kinematics-dependent couplings.
540 0 : int idAbs1 = abs(id1);
541 0 : int idAbs2 = abs(id2);
542 0 : int i3 = abs(id3chi);
543 0 : int i4 = abs(id4chi);
544 :
545 : // Flavour-dependent kinematics-dependent couplings.
546 0 : complex QuLL(0.0),QtLL(0.0),QuRR(0.0),QtRR(0.0);
547 0 : complex QuLR(0.0),QtLR(0.0),QuRL(0.0),QtRL(0.0);
548 :
549 : double *LqqZloc;
550 : double *RqqZloc;
551 : complex (*LsduXloc)[4][3];
552 : complex (*RsduXloc)[4][3];
553 : complex (*LsudXloc)[4][3];
554 : complex (*RsudXloc)[4][3];
555 :
556 : int iShift(0);
557 0 : if( idAbs1 > 10 && idAbs1 < 17 ) {
558 : iShift+=10;
559 0 : LqqZloc = coupSUSYPtr->LllZ;
560 0 : RqqZloc = coupSUSYPtr->RllZ;
561 0 : LsduXloc = coupSUSYPtr->LslvX;
562 0 : RsduXloc = coupSUSYPtr->RslvX;
563 0 : LsudXloc = coupSUSYPtr->LsvlX;
564 0 : RsudXloc = coupSUSYPtr->RsvlX;
565 0 : } else {
566 0 : LqqZloc = coupSUSYPtr->LqqZ;
567 0 : RqqZloc = coupSUSYPtr->RqqZ;
568 0 : LsduXloc = coupSUSYPtr->LsduX;
569 0 : RsduXloc = coupSUSYPtr->RsduX;
570 0 : LsudXloc = coupSUSYPtr->LsudX;
571 0 : RsudXloc = coupSUSYPtr->RsudX;
572 : }
573 :
574 : // Add Z/gamma* for same-flavour in-quarks
575 0 : if (idAbs1 == idAbs2) {
576 :
577 0 : QuLL = -LqqZloc[idAbs1-iShift]*conj(coupSUSYPtr->ORp[i3][i4]);
578 0 : QtLL = -LqqZloc[idAbs1-iShift]*conj(coupSUSYPtr->OLp[i3][i4]);
579 0 : QuRR = -RqqZloc[idAbs1-iShift]*conj(coupSUSYPtr->OLp[i3][i4]);
580 0 : QtRR = -RqqZloc[idAbs1-iShift]*conj(coupSUSYPtr->ORp[i3][i4]);
581 :
582 0 : QuLL *= propZ / 2.0 / (1.0-coupSUSYPtr->sin2W);
583 0 : QtLL *= propZ / 2.0 / (1.0-coupSUSYPtr->sin2W);
584 0 : QuRR *= propZ / 2.0 / (1.0-coupSUSYPtr->sin2W);
585 0 : QtRR *= propZ / 2.0 / (1.0-coupSUSYPtr->sin2W);
586 :
587 : // s-channel gamma* (only for same-type charginos)
588 0 : if (i3 == i4) {
589 :
590 : // Charge of in-particles
591 0 : double q = particleDataPtr->chargeType(idAbs1)/3.0;
592 0 : QuLL += q * coupSUSYPtr->sin2W / sH;
593 0 : QuRR += q * coupSUSYPtr->sin2W / sH;
594 0 : QtLL += q * coupSUSYPtr->sin2W / sH;
595 0 : QtRR += q * coupSUSYPtr->sin2W / sH;
596 :
597 0 : }
598 : }
599 :
600 0 : int iG1 = (abs(id1)+1-iShift)/2;
601 0 : int iG2 = (abs(id2)+1-iShift)/2;
602 :
603 : // Add t- or u-channel squark flavour sums to QmXY couplings
604 0 : for (int ksq=1; ksq<=6; ksq++) {
605 :
606 0 : if(id1 % 2 == 0) {
607 :
608 : // u-channel diagrams only
609 : // up-type incoming; u-channel ~d
610 :
611 0 : int idsd = ((ksq+2)/3)*1000000 + 2*((ksq-1) % 3) + 1;
612 0 : idsd +=iShift;
613 0 : double msq = particleDataPtr->m0(idsd);
614 0 : double ufac = 2.0 * (uH - pow2(msq));
615 :
616 : //u-ubar -> chi-chi+
617 0 : QuLL += LsduXloc[ksq][iG2][i3]
618 0 : * conj(LsduXloc[ksq][iG1][i4]) / ufac;
619 0 : QuRR += RsduXloc[ksq][iG2][i3]
620 0 : * conj(RsduXloc[ksq][iG1][i4]) / ufac;
621 0 : QuLR += RsduXloc[ksq][iG2][i3]
622 0 : * conj(LsduXloc[ksq][iG1][i4]) / ufac;
623 0 : QuRL += LsduXloc[ksq][iG2][i3]
624 0 : * conj(RsduXloc[ksq][iG1][i4]) / ufac;
625 :
626 0 : }else{
627 : // t-channel diagrams only;
628 : // down-type incoming; t-channel ~u
629 :
630 0 : int idsu = ((ksq+2)/3)*1000000 + 2*((ksq-1) % 3) + 2;
631 0 : idsu += iShift;
632 0 : double msq = particleDataPtr->m0(idsu);
633 0 : double tfac = 2.0 * (tH - pow2(msq));
634 :
635 : //d-dbar -> chi-chi+
636 0 : QtLL -= LsudXloc[ksq][iG1][i3]
637 0 : * conj(LsudXloc[ksq][iG2][i4]) / tfac;
638 0 : QtRR -= RsudXloc[ksq][iG1][i3]
639 0 : * conj(RsudXloc[ksq][iG2][i4]) / tfac;
640 0 : QtLR += LsudXloc[ksq][iG1][i3]
641 0 : * conj(RsudXloc[ksq][iG2][i4]) / tfac;
642 0 : QtRL += RsudXloc[ksq][iG1][i3]
643 0 : * conj(LsudXloc[ksq][iG2][i4]) / tfac;
644 :
645 0 : }
646 : }
647 : // Compute matrix element weight
648 : double weight = 0;
649 :
650 : // Average over separate helicity contributions
651 : // LL (ha = -1, hb = +1) (divided by 4 for average)
652 0 : weight += norm(QuLL) * ui * uj + norm(QtLL) * ti * tj
653 0 : + 2 * real(conj(QuLL) * QtLL) * m3 * m4 * sH;
654 : // RR (ha = 1, hb = -1) (divided by 4 for average)
655 0 : weight += norm(QtRR) * ti * tj + norm(QuRR) * ui * uj
656 0 : + 2 * real(conj(QuRR) * QtRR) * m3 * m4 * sH;
657 : // RL (ha = 1, hb = 1) (divided by 4 for average)
658 0 : weight += norm(QuRL) * ui * uj + norm(QtRL) * ti * tj
659 0 : + real(conj(QuRL) * QtRL) * (uH * tH - s3 * s4);
660 : // LR (ha = -1, hb = -1) (divided by 4 for average)
661 0 : weight += norm(QuLR) * ui * uj + norm(QtLR) * ti * tj
662 0 : + real(conj(QuLR) * QtLR) * (uH * tH - s3 * s4);
663 :
664 0 : double colorFactor = ( idAbs1 > 10 && idAbs1 < 17 ) ? 3.0 : 1.0;
665 :
666 : // Cross section, including colour factor.
667 0 : double sigma = sigma0 * weight * colorFactor;
668 :
669 : // Answer.
670 : return sigma;
671 :
672 0 : }
673 :
674 : //==========================================================================
675 :
676 : // Sigma2qgchi0squark
677 : // Cross section for gaugino-squark production: neutralino-squark
678 :
679 : //--------------------------------------------------------------------------
680 :
681 : // Initialize process.
682 :
683 : void Sigma2qg2chi0squark::initProc() {
684 :
685 : //Typecast to the correct couplings
686 0 : coupSUSYPtr = (CoupSUSY*) couplingsPtr;
687 :
688 : // Construct name of process.
689 0 : if (id4 % 2 == 0) {
690 0 : nameSave = "q g -> " + particleDataPtr->name(id3) + " "
691 0 : + particleDataPtr->name(id4) + " + c.c. (q=u,c)";
692 0 : }
693 : else {
694 0 : nameSave = "q g -> " + particleDataPtr->name(id3) + " "
695 0 : + particleDataPtr->name(id4) + " + c.c. (q=d,s,b)";
696 : }
697 :
698 : // Secondary open width fraction.
699 0 : openFracPair = particleDataPtr->resOpenFrac(id3, id4);
700 :
701 0 : }
702 :
703 : //--------------------------------------------------------------------------
704 :
705 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
706 :
707 : void Sigma2qg2chi0squark::sigmaKin() {
708 :
709 : // Common flavour-independent factor.
710 : // tmp: alphaS = 0.1 for counter-checks
711 0 : double nChi = 6.0 * coupSUSYPtr->sin2W * (1 - coupSUSYPtr->sin2W);
712 0 : sigma0 = M_PI / sH2 / nChi * alpEM * alpS
713 0 : * openFracPair;
714 :
715 : // Auxiliary factors for use below
716 0 : ui = uH - s3;
717 0 : uj = uH - s4;
718 0 : ti = tH - s3;
719 0 : tj = tH - s4;
720 :
721 0 : }
722 :
723 : //--------------------------------------------------------------------------
724 :
725 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
726 :
727 : double Sigma2qg2chi0squark::sigmaHat() {
728 :
729 : // Antiquark -> antisquark
730 0 : int idq = id1;
731 0 : if (id1 == 21 || id1 == 22) idq = id2;
732 0 : if (idq < 0) {
733 0 : id4 = -abs(id4);
734 0 : } else {
735 0 : id4 = abs(id4);
736 : }
737 :
738 : // tmp: only allow incoming quarks on side 1
739 : // if (id1 < 0 || id1 == 21) return 0.0;
740 :
741 : // Generation index
742 0 : int iGq = (abs(idq)+1)/2;
743 :
744 : // Only accept u(bar) -> ~u(bar) and d(bar) -> ~d(bar)
745 0 : if (particleDataPtr->chargeType(idq) != particleDataPtr->chargeType(id4))
746 0 : return 0.0;
747 :
748 : // Couplings
749 0 : complex LsqqX, RsqqX;
750 0 : if (idq % 2 == 0) {
751 0 : LsqqX = coupSUSYPtr->LsuuX[id4sq][iGq][id3chi];
752 0 : RsqqX = coupSUSYPtr->RsuuX[id4sq][iGq][id3chi];
753 0 : }
754 : else {
755 0 : LsqqX = coupSUSYPtr->LsddX[id4sq][iGq][id3chi];
756 0 : RsqqX = coupSUSYPtr->RsddX[id4sq][iGq][id3chi];
757 : }
758 :
759 : // Prefactors : swap u and t if gq instead of qg
760 : double fac1, fac2;
761 0 : if (idq == id1) {
762 0 : fac1 = -ui/sH + 2.0 * ( uH*tH - s4*s3 )/sH/tj;
763 0 : fac2 = ti/tj * ( (tH + s4)/tj + (ti - uj)/sH );
764 0 : } else {
765 0 : fac1 = -ti/sH + 2.0 * ( uH*tH - s4*s3 )/sH/uj;
766 0 : fac2 = ui/uj * ( (uH + s4)/uj + (ui - tj)/sH );
767 : }
768 :
769 : // Compute matrix element weight
770 : double weight = 0.0;
771 :
772 : // Average over separate helicity contributions
773 : // (for qbar g : ha -> -ha )
774 : // LL (ha = -1, hb = +1) (divided by 4 for average)
775 0 : weight += fac2 * norm(LsqqX) / 2.0;
776 : // RR (ha = 1, hb = -1) (divided by 4 for average)
777 0 : weight += fac2 * norm(RsqqX) / 2.0;
778 : // RL (ha = 1, hb = 1) (divided by 4 for average)
779 0 : weight += fac2 * norm(RsqqX) / 2.0 + fac1 * norm(RsqqX);
780 : // LR (ha = -1, hb = -1) (divided by 4 for average)
781 0 : weight += fac2 * norm(LsqqX) / 2.0 + fac1 * norm(LsqqX);
782 :
783 0 : double sigma = sigma0 * weight;
784 :
785 : // Answer.
786 : return sigma;
787 :
788 0 : }
789 :
790 : //--------------------------------------------------------------------------
791 :
792 : // Select identity, colour and anticolour.
793 :
794 : void Sigma2qg2chi0squark::setIdColAcol() {
795 :
796 : // Set flavours.
797 0 : setId( id1, id2, id3, (id1*id2 > 0 ? abs(id4) : -abs(id4)));
798 :
799 : // Colour flow topology. Swap if first is gluon, or when antiquark.
800 0 : if (id1 != 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
801 0 : else setColAcol( 1, 2, 2, 0, 0, 0, 1, 0);
802 0 : if (id1*id2 < 0) swapColAcol();
803 :
804 0 : }
805 :
806 : //==========================================================================
807 :
808 : // Sigma2qg2charsquark
809 : // Cross section for gaugino-squark production: chargino-squark
810 :
811 : //--------------------------------------------------------------------------
812 :
813 : // Initialize process.
814 :
815 : void Sigma2qg2charsquark::initProc() {
816 :
817 : //Typecast to the correct couplings
818 0 : coupSUSYPtr = (CoupSUSY*) couplingsPtr;
819 :
820 : // Construct name of process.
821 0 : if (id4 % 2 == 0) {
822 0 : nameSave = "q g -> " + particleDataPtr->name(id3) + " "
823 0 : + particleDataPtr->name(id4) + " + c.c. (q=d,s,b)";
824 0 : }
825 : else {
826 0 : nameSave = "q g -> " + particleDataPtr->name(id3) + " "
827 0 : + particleDataPtr->name(id4) + " + c.c. (q=u,c)";
828 : }
829 :
830 : // Secondary open width fraction.
831 0 : openFracPair = particleDataPtr->resOpenFrac(id3Sav, id4Sav);
832 :
833 0 : }
834 :
835 : //--------------------------------------------------------------------------
836 :
837 : void Sigma2qg2charsquark::sigmaKin() {
838 :
839 : // Common flavour-independent factor.
840 : // tmp: alphaS = 0.1 for counter-checks
841 0 : double nChi = 12.0 * coupSUSYPtr->sin2W;
842 0 : sigma0 = M_PI / sH2 / nChi * alpEM * alpS
843 0 : * openFracPair;
844 :
845 : // Auxiliary factors for use below
846 0 : ui = uH - s3;
847 0 : uj = uH - s4;
848 0 : ti = tH - s3;
849 0 : tj = tH - s4;
850 :
851 0 : }
852 :
853 : //--------------------------------------------------------------------------
854 :
855 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
856 :
857 : double Sigma2qg2charsquark::sigmaHat() {
858 :
859 : // Antiquark -> antisquark
860 0 : int idq = (id1 == 21) ? id2 : id1;
861 0 : if (idq > 0) {
862 0 : id3 = id3Sav;
863 0 : id4 = id4Sav;
864 0 : } else {
865 0 : id3 = -id3Sav;
866 0 : id4 = -id4Sav;
867 : }
868 :
869 : // Only accept u(bar) -> ~d(bar) and d(bar) -> ~u(bar)
870 0 : if (particleDataPtr->chargeType(idq) == particleDataPtr->chargeType(id4))
871 0 : return 0.0;
872 :
873 : // Generation index
874 0 : int iGq = (abs(idq)+1)/2;
875 :
876 : // Couplings
877 0 : complex LsqqX, RsqqX;
878 0 : if (idq % 2 == 0) {
879 0 : LsqqX = coupSUSYPtr->LsduX[id4sq][iGq][id3chi];
880 0 : RsqqX = coupSUSYPtr->RsduX[id4sq][iGq][id3chi];
881 0 : }
882 : else {
883 0 : LsqqX = coupSUSYPtr->LsudX[id4sq][iGq][id3chi];
884 0 : RsqqX = coupSUSYPtr->RsudX[id4sq][iGq][id3chi];
885 : }
886 :
887 : // Prefactors : swap u and t if gq instead of qg
888 : double fac1, fac2;
889 0 : if (idq == id1) {
890 0 : fac1 = -ui/sH + 2.0 * ( uH*tH - s4*s3 )/sH/tj;
891 0 : fac2 = ti/tj * ( (tH + s4)/tj + (ti - uj)/sH );
892 0 : } else {
893 0 : fac1 = -ti/sH + 2.0 * ( uH*tH - s4*s3 )/sH/uj;
894 0 : fac2 = ui/uj * ( (uH + s4)/uj + (ui - tj)/sH );
895 : }
896 :
897 : // Compute matrix element weight
898 : double weight = 0.0;
899 :
900 : // Average over separate helicity contributions
901 : // (a, b refers to qg configuration)
902 : // LL (ha = -1, hb = +1) (divided by 4 for average)
903 0 : weight += fac2 * norm(LsqqX) / 2.0;
904 : // RR (ha = 1, hb = -1) (divided by 4 for average)
905 0 : weight += fac2 * norm(RsqqX) / 2.0;
906 : // RL (ha = 1, hb = 1) (divided by 4 for average)
907 0 : weight += fac2 * norm(RsqqX) / 2.0 + fac1 * norm(RsqqX);
908 : // LR (ha = -1, hb = -1) (divided by 4 for average)
909 0 : weight += fac2 * norm(LsqqX) / 2.0 + fac1 * norm(LsqqX);
910 :
911 0 : double sigma = sigma0 * weight;
912 :
913 : // Answer.
914 0 : return sigma * openFracPair;
915 :
916 0 : }
917 :
918 : //--------------------------------------------------------------------------
919 :
920 : // Select identity, colour and anticolour.
921 :
922 : void Sigma2qg2charsquark::setIdColAcol() {
923 :
924 : // Set flavours.
925 0 : if (id1 > 0 && id2 > 0) {
926 0 : setId( id1, id2, id3Sav, id4Sav);
927 0 : } else {
928 0 : setId( id1, id2,-id3Sav,-id4Sav);
929 : }
930 :
931 : // Colour flow topology. Swap if first is gluon, or when antiquark.
932 0 : if (id1 != 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
933 0 : else setColAcol( 1, 2, 2, 0, 0, 0, 1, 0);
934 0 : if (id1 < 0 || id2 < 0) swapColAcol();
935 :
936 0 : }
937 :
938 : //==========================================================================
939 :
940 : // Sigma2qq2squarksquark
941 : // Cross section for squark-squark production
942 :
943 : //--------------------------------------------------------------------------
944 :
945 : // Initialize process.
946 :
947 : void Sigma2qq2squarksquark::initProc() {
948 :
949 : //Typecast to the correct couplings
950 0 : coupSUSYPtr = (CoupSUSY*) couplingsPtr;
951 :
952 : // Extract mass-ordering indices
953 0 : iGen3 = 3*(abs(id3Sav)/2000000) + (abs(id3Sav)%10+1)/2;
954 0 : iGen4 = 3*(abs(id4Sav)/2000000) + (abs(id4Sav)%10+1)/2;
955 :
956 : // Is this a ~u_i ~d_j fial state or ~d_i ~d_j, ~u_i ~u_j
957 0 : if (abs(id3Sav) % 2 == abs(id4Sav) % 2) isUD = false;
958 0 : else isUD = true;
959 :
960 : // Derive name
961 0 : nameSave = "q q' -> "+particleDataPtr->name(abs(id3Sav))+" "
962 0 : +particleDataPtr->name(abs(id4Sav))+" + c.c.";
963 :
964 : // Count 5 neutralinos in NMSSM
965 0 : nNeut = (coupSUSYPtr->isNMSSM ? 5 : 4);
966 :
967 : // Store mass squares of all possible internal propagator lines
968 0 : m2Glu = pow2(particleDataPtr->m0(1000021));
969 0 : m2Neut.resize(nNeut+1);
970 0 : for (int iNeut=1;iNeut<=nNeut;iNeut++) {
971 0 : m2Neut[iNeut] = pow2(particleDataPtr->m0(coupSUSYPtr->idNeut(iNeut)));
972 : }
973 0 : m2Char.resize(3);
974 0 : m2Char[1] = pow2(particleDataPtr->m0(coupSUSYPtr->idChar(1)));
975 0 : m2Char[2] = pow2(particleDataPtr->m0(coupSUSYPtr->idChar(2)));
976 :
977 : // Set sizes of some arrays to be used below
978 0 : tNeut.resize(nNeut+1);
979 0 : uNeut.resize(nNeut+1);
980 0 : tChar.resize(3);
981 0 : uChar.resize(3);
982 :
983 : // Secondary open width fraction.
984 0 : openFracPair = particleDataPtr->resOpenFrac(id3Sav, id4Sav);
985 :
986 : // Selection of interference terms
987 0 : onlyQCD = settingsPtr->flag("SUSY:qq2squarksquark:onlyQCD");
988 0 : }
989 :
990 : //--------------------------------------------------------------------------
991 :
992 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
993 :
994 : void Sigma2qq2squarksquark::sigmaKin() {
995 :
996 : // Weak mixing
997 0 : double xW = coupSUSYPtr->sin2W;
998 :
999 : // pi/sH2
1000 0 : double comFacHat = M_PI/sH2 * openFracPair;
1001 :
1002 : // Channel-dependent but flavor-independent pre-factors
1003 0 : sigmaNeut = comFacHat * pow2(alpEM) / pow2(xW) / pow2(1-xW);
1004 0 : sigmaGlu = comFacHat * 2.0 * pow2(alpS) / 9.0;
1005 0 : if (isUD) {
1006 0 : sigmaChar = comFacHat * pow2(alpEM) / 4.0 / pow2(xW);
1007 0 : sigmaCharNeut = comFacHat * pow2(alpEM) / 3.0 / pow2(xW) / (1-xW);
1008 0 : sigmaCharGlu = comFacHat * 4.0 * alpEM * alpS / 9.0 / xW;
1009 0 : sigmaNeutGlu = 0.0;
1010 0 : } else {
1011 0 : sigmaChar = 0.0;
1012 0 : sigmaCharNeut = 0.0;
1013 0 : sigmaCharGlu = 0.0;
1014 0 : sigmaNeutGlu = comFacHat * 8.0 * alpEM * alpS / 9.0 / xW/(1-xW);
1015 : }
1016 :
1017 0 : }
1018 :
1019 : //--------------------------------------------------------------------------
1020 :
1021 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1022 :
1023 : double Sigma2qq2squarksquark::sigmaHat() {
1024 :
1025 : // In-pair must be same-sign
1026 0 : if (id1 * id2 < 0) return 0.0;
1027 :
1028 : // Check correct charge sum
1029 0 : if (isUD && abs(id1) %2 == abs(id2) % 2) return 0.0;
1030 0 : if (!isUD && abs(id1) % 2 != abs(id2) % 2) return 0.0;
1031 0 : if (!isUD && abs(id1) % 2 != abs(id3Sav) % 2) return 0.0;
1032 :
1033 : // Coded sigma is for ud -> ~q~q'. Swap t and u for du -> ~q~q'.
1034 0 : swapTU = (isUD && abs(id1) % 2 == 0);
1035 0 : int idIn1A = (swapTU) ? abs(id2) : abs(id1);
1036 0 : int idIn2A = (swapTU) ? abs(id1) : abs(id2);
1037 :
1038 : // Auxiliary factors for use below
1039 0 : tGlu = tH - m2Glu;
1040 0 : uGlu = uH - m2Glu;
1041 0 : for (int i=1; i<= nNeut; i++) {
1042 0 : tNeut[i] = tH - m2Neut[i];
1043 0 : uNeut[i] = uH - m2Neut[i];
1044 0 : if (isUD && i <= 2) {
1045 0 : tChar[i] = tH - m2Char[i];
1046 0 : uChar[i] = uH - m2Char[i];
1047 0 : }
1048 : }
1049 :
1050 : // Generation indices of incoming particles
1051 0 : int iGen1 = (abs(idIn1A)+1)/2;
1052 0 : int iGen2 = (abs(idIn2A)+1)/2;
1053 :
1054 : // Initial values for pieces used for color-flow selection below
1055 0 : sumCt = 0.0;
1056 0 : sumCu = 0.0;
1057 0 : sumNt = 0.0;
1058 0 : sumNu = 0.0;
1059 0 : sumGt = 0.0;
1060 0 : sumGu = 0.0;
1061 0 : sumInterference = 0.0;
1062 :
1063 : // Common factor for LR and RL contributions
1064 0 : double facTU = uH*tH-s3*s4;
1065 :
1066 : // Case A) Opposite-isospin: qq' -> ~d~u
1067 0 : if ( isUD ) {
1068 :
1069 : // t-channel Charginos
1070 0 : for (int k=1;k<=2;k++) {
1071 :
1072 : // Skip if only including gluinos
1073 0 : if (onlyQCD) continue;
1074 :
1075 0 : for (int l=1;l<=2;l++) {
1076 :
1077 : // kl-dependent factor for LL and RR contributions
1078 0 : double facMS = sH*sqrt(m2Char[k]*m2Char[l]);
1079 :
1080 : // Note: Ckl defined as in [Boz07] with sigmaChar factored out
1081 : // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
1082 0 : complex Ckl[3][3];
1083 0 : Ckl[1][1] = facMS
1084 0 : * coupSUSYPtr->LsudX[iGen4][iGen2][k]
1085 0 : * conj(coupSUSYPtr->LsduX[iGen3][iGen1][k])
1086 0 : * conj(coupSUSYPtr->LsudX[iGen4][iGen2][l])
1087 0 : * coupSUSYPtr->LsduX[iGen3][iGen1][l];
1088 0 : Ckl[1][2] = facTU
1089 0 : * coupSUSYPtr->RsudX[iGen4][iGen2][k]
1090 0 : * conj(coupSUSYPtr->LsduX[iGen3][iGen1][k])
1091 0 : * conj(coupSUSYPtr->RsudX[iGen4][iGen2][l])
1092 0 : * coupSUSYPtr->LsduX[iGen3][iGen1][l];
1093 0 : Ckl[2][1] = facTU
1094 0 : * coupSUSYPtr->LsudX[iGen4][iGen2][k]
1095 0 : * conj(coupSUSYPtr->RsduX[iGen3][iGen1][k])
1096 0 : * conj(coupSUSYPtr->LsudX[iGen4][iGen2][l])
1097 0 : * coupSUSYPtr->RsduX[iGen3][iGen1][l];
1098 0 : Ckl[2][2] = facMS
1099 0 : * coupSUSYPtr->RsudX[iGen4][iGen2][k]
1100 0 : * conj(coupSUSYPtr->RsduX[iGen3][iGen1][k])
1101 0 : * conj(coupSUSYPtr->RsudX[iGen4][iGen2][l])
1102 0 : * coupSUSYPtr->RsduX[iGen3][iGen1][l];
1103 :
1104 : // Add to sum of t-channel charginos
1105 0 : sumCt += sigmaChar * real(Ckl[1][1] + Ckl[1][2] + Ckl[2][1]
1106 0 : + Ckl[2][2]) / tChar[k] / tChar[l];
1107 :
1108 0 : }
1109 0 : }
1110 :
1111 : // Skip if only including gluinos
1112 0 : if (!onlyQCD) {
1113 :
1114 : // u-channel Neutralinos
1115 0 : for (int k=1;k<=nNeut;k++) {
1116 0 : for (int l=1;l<=nNeut;l++) {
1117 :
1118 : // kl-dependent factor for LL, RR contributions
1119 0 : double facMS = sH*sqrt(m2Neut[k]*m2Neut[l]);
1120 :
1121 : // Note: Nkl defined as in [Boz07] with sigmaNeut factored out
1122 : // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
1123 0 : complex Nkl[3][3];
1124 0 : Nkl[1][1] = facMS
1125 0 : * conj(coupSUSYPtr->LsuuX[iGen4][iGen1][k])
1126 0 : * conj(coupSUSYPtr->LsddX[iGen3][iGen2][k])
1127 0 : * coupSUSYPtr->LsuuX[iGen4][iGen1][l]
1128 0 : * coupSUSYPtr->LsddX[iGen3][iGen2][l];
1129 0 : Nkl[1][2] = facTU
1130 0 : * conj(coupSUSYPtr->LsuuX[iGen4][iGen1][k])
1131 0 : * conj(coupSUSYPtr->RsddX[iGen3][iGen2][k])
1132 0 : * coupSUSYPtr->LsuuX[iGen4][iGen1][l]
1133 0 : * coupSUSYPtr->RsddX[iGen3][iGen2][l];
1134 0 : Nkl[2][1] = facTU
1135 0 : * conj(coupSUSYPtr->RsuuX[iGen4][iGen1][k])
1136 0 : * conj(coupSUSYPtr->LsddX[iGen3][iGen2][k])
1137 0 : * coupSUSYPtr->RsuuX[iGen4][iGen1][l]
1138 0 : * coupSUSYPtr->LsddX[iGen3][iGen2][l];
1139 0 : Nkl[2][2] = facMS
1140 0 : * conj(coupSUSYPtr->RsuuX[iGen4][iGen1][k])
1141 0 : * conj(coupSUSYPtr->RsddX[iGen3][iGen2][k])
1142 0 : * coupSUSYPtr->RsuuX[iGen4][iGen1][l]
1143 0 : * coupSUSYPtr->RsddX[iGen3][iGen2][l];
1144 :
1145 : // Add to sum of u-channel neutralinos
1146 0 : sumNu += sigmaNeut / uNeut[k] / uNeut[l]
1147 0 : * real(Nkl[1][1] + Nkl[1][2] + Nkl[2][1] + Nkl[2][2]);
1148 :
1149 0 : }
1150 : }
1151 0 : }
1152 :
1153 : // u-channel gluino
1154 : // Note: Gij defined as in [Boz07] with sigmaGlu factored out
1155 : // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
1156 : double Gij[3][3];
1157 0 : Gij[1][1] = norm(coupSUSYPtr->LsuuG[iGen4][iGen1]
1158 0 : * coupSUSYPtr->LsddG[iGen3][iGen2]);
1159 0 : Gij[1][2] = norm(coupSUSYPtr->LsuuG[iGen4][iGen1]
1160 0 : * coupSUSYPtr->RsddG[iGen3][iGen2]);
1161 0 : Gij[2][1] = norm(coupSUSYPtr->RsuuG[iGen4][iGen1]
1162 0 : * coupSUSYPtr->LsddG[iGen3][iGen2]);
1163 0 : Gij[2][2] = norm(coupSUSYPtr->RsuuG[iGen4][iGen1]
1164 0 : * coupSUSYPtr->RsddG[iGen3][iGen2]);
1165 0 : Gij[1][1] *= sH*m2Glu;
1166 0 : Gij[1][2] *= facTU;
1167 0 : Gij[2][1] *= facTU;
1168 0 : Gij[2][2] *= sH*m2Glu;
1169 : // Sum over polarizations
1170 0 : sumGu += sigmaGlu * (Gij[1][1] + Gij[1][2] + Gij[2][1] + Gij[2][2])
1171 0 : / pow2(uGlu);
1172 :
1173 :
1174 : // EW Interference terms: Skip if only including gluinos
1175 0 : if (!onlyQCD) {
1176 :
1177 : // chargino-neutralino interference
1178 0 : for (int k=1;k<=2;k++) {
1179 0 : for (int l=1;l<=nNeut;l++) {
1180 : // Note: CNkl defined as in [Boz07] with pi/sH2 factored out
1181 : // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
1182 : double CNkl[3][3];
1183 0 : CNkl[1][1] = real(coupSUSYPtr->LsudX[iGen4][iGen2][k]
1184 0 : * conj(coupSUSYPtr->LsduX[iGen3][iGen1][k])
1185 0 : * coupSUSYPtr->LsuuX[iGen4][iGen1][l]
1186 0 : * coupSUSYPtr->LsddX[iGen3][iGen2][l]);
1187 0 : CNkl[1][2] = real(coupSUSYPtr->RsudX[iGen4][iGen2][k]
1188 0 : * conj(coupSUSYPtr->LsduX[iGen3][iGen1][k])
1189 0 : * coupSUSYPtr->LsuuX[iGen4][iGen1][l]
1190 0 : * coupSUSYPtr->RsddX[iGen3][iGen2][l]);
1191 0 : CNkl[2][1] = real(coupSUSYPtr->LsudX[iGen4][iGen2][k]
1192 0 : * conj(coupSUSYPtr->RsduX[iGen3][iGen1][k])
1193 0 : * coupSUSYPtr->RsuuX[iGen4][iGen1][l]
1194 0 : * coupSUSYPtr->LsddX[iGen3][iGen2][l]);
1195 0 : CNkl[2][2] = real(coupSUSYPtr->RsudX[iGen4][iGen2][k]
1196 0 : * conj(coupSUSYPtr->RsduX[iGen3][iGen1][k])
1197 0 : * coupSUSYPtr->RsuuX[iGen4][iGen1][l]
1198 0 : * coupSUSYPtr->RsddX[iGen3][iGen2][l]);
1199 0 : CNkl[1][1] *= sH*sqrt(m2Char[k]*m2Neut[l]);
1200 0 : CNkl[1][2] *= uH*tH-s3*s4;
1201 0 : CNkl[2][1] *= uH*tH-s3*s4;
1202 0 : CNkl[2][2] *= sH*sqrt(m2Char[k]*m2Neut[l]);
1203 : // Sum over polarizations
1204 0 : sumInterference += sigmaCharNeut * (CNkl[1][1] + CNkl[1][2]
1205 0 : + CNkl[2][1] + CNkl[2][2]) / tChar[k] / uNeut[l];
1206 : }
1207 : }
1208 :
1209 : // chargino-gluino interference
1210 0 : for (int k=1;k<=2;k++) {
1211 : // Note: CGk defined as in [Boz07] with sigmaCharGlu factored out
1212 : // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
1213 : double CGk[3][3];
1214 0 : CGk[1][1] = real(coupSUSYPtr->LsudX[iGen4][iGen2][k]
1215 0 : * conj(coupSUSYPtr->LsduX[iGen3][iGen1][k])
1216 0 : * conj(coupSUSYPtr->LsuuG[iGen4][iGen1])
1217 0 : * conj(coupSUSYPtr->LsddG[iGen3][iGen2]));
1218 0 : CGk[1][2] = real(coupSUSYPtr->RsudX[iGen4][iGen2][k]
1219 0 : * conj(coupSUSYPtr->LsduX[iGen3][iGen1][k])
1220 0 : * conj(coupSUSYPtr->LsuuG[iGen4][iGen1])
1221 0 : * conj(coupSUSYPtr->RsddG[iGen3][iGen2]));
1222 0 : CGk[2][1] = real(coupSUSYPtr->LsudX[iGen4][iGen2][k]
1223 0 : * conj(coupSUSYPtr->RsduX[iGen3][iGen1][k])
1224 0 : * conj(coupSUSYPtr->RsuuG[iGen4][iGen1])
1225 0 : * conj(coupSUSYPtr->LsddG[iGen3][iGen2]));
1226 0 : CGk[2][2] = real(coupSUSYPtr->RsudX[iGen4][iGen2][k]
1227 0 : * conj(coupSUSYPtr->RsduX[iGen3][iGen1][k])
1228 0 : * conj(coupSUSYPtr->RsuuG[iGen4][iGen1])
1229 0 : * conj(coupSUSYPtr->RsddG[iGen3][iGen2]));
1230 0 : CGk[1][1] *= sH*sqrt(m2Glu*m2Char[k]);
1231 0 : CGk[1][2] *= uH*tH-s3*s4;
1232 0 : CGk[2][1] *= uH*tH-s3*s4;
1233 0 : CGk[2][2] *= sH*sqrt(m2Glu*m2Char[k]);
1234 : // Sum over polarizations
1235 0 : sumInterference += sigmaGlu * (CGk[1][1] + CGk[1][2] + CGk[2][1]
1236 0 : + CGk[2][2]) / uGlu / tChar[k];
1237 : }
1238 0 : }
1239 0 : }
1240 :
1241 : // Case B) Same-isospin: qq' -> ~d~d , ~u~u
1242 : else {
1243 :
1244 : // t-channel + u-channel Neutralinos + t/u interference
1245 : // Skip if only including gluinos
1246 0 : if (!onlyQCD) {
1247 0 : for (int k=1;k<=nNeut;k++) {
1248 0 : for (int l=1;l<=nNeut;l++) {
1249 :
1250 : // kl-dependent factor for LL and RR contributions
1251 0 : double facMS = sH * particleDataPtr->m0(coupSUSYPtr->idNeut(k))
1252 0 : * particleDataPtr->m0(coupSUSYPtr->idNeut(l));
1253 :
1254 : // Note: Nxkl defined as in [Boz07] with sigmaNeut factored out
1255 : // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
1256 0 : complex NTkl[3][3], NUkl[3][3], NTUkl[3][3];
1257 0 : NTkl[1][1] = facMS
1258 0 : * conj(coupSUSYPtr->getLsqqX(iGen4,idIn2A,k))
1259 0 : * conj(coupSUSYPtr->getLsqqX(iGen3,idIn1A,k))
1260 0 : * coupSUSYPtr->getLsqqX(iGen4,idIn2A,l)
1261 0 : * coupSUSYPtr->getLsqqX(iGen3,idIn1A,l);
1262 0 : NTkl[1][2] = facTU
1263 0 : * conj(coupSUSYPtr->getRsqqX(iGen4,idIn2A,k))
1264 0 : * conj(coupSUSYPtr->getLsqqX(iGen3,idIn1A,k))
1265 0 : * coupSUSYPtr->getRsqqX(iGen4,idIn2A,l)
1266 0 : * coupSUSYPtr->getLsqqX(iGen3,idIn1A,l);
1267 0 : NTkl[2][1] = facTU
1268 0 : * conj(coupSUSYPtr->getLsqqX(iGen4,idIn2A,k))
1269 0 : * conj(coupSUSYPtr->getRsqqX(iGen3,idIn1A,k))
1270 0 : * coupSUSYPtr->getLsqqX(iGen4,idIn2A,l)
1271 0 : * coupSUSYPtr->getRsqqX(iGen3,idIn1A,l);
1272 0 : NTkl[2][2] = facMS
1273 0 : * conj(coupSUSYPtr->getRsqqX(iGen4,idIn2A,k))
1274 0 : * conj(coupSUSYPtr->getRsqqX(iGen3,idIn1A,k))
1275 0 : * coupSUSYPtr->getRsqqX(iGen4,idIn2A,l)
1276 0 : * coupSUSYPtr->getRsqqX(iGen3,idIn1A,l);
1277 0 : NUkl[1][1] = facMS
1278 0 : * conj(coupSUSYPtr->getLsqqX(iGen3,idIn2A,k))
1279 0 : * conj(coupSUSYPtr->getLsqqX(iGen4,idIn1A,k))
1280 0 : * coupSUSYPtr->getLsqqX(iGen3,idIn2A,l)
1281 0 : * coupSUSYPtr->getLsqqX(iGen4,idIn1A,l);
1282 0 : NUkl[1][2] = facTU
1283 0 : * conj(coupSUSYPtr->getRsqqX(iGen3,idIn2A,k))
1284 0 : * conj(coupSUSYPtr->getLsqqX(iGen4,idIn1A,k))
1285 0 : * coupSUSYPtr->getRsqqX(iGen3,idIn2A,l)
1286 0 : * coupSUSYPtr->getLsqqX(iGen4,idIn1A,l);
1287 0 : NUkl[2][1] = facTU
1288 0 : * conj(coupSUSYPtr->getLsqqX(iGen3,idIn2A,k))
1289 0 : * conj(coupSUSYPtr->getRsqqX(iGen4,idIn1A,k))
1290 0 : * coupSUSYPtr->getLsqqX(iGen3,idIn2A,l)
1291 0 : * coupSUSYPtr->getRsqqX(iGen4,idIn1A,l);
1292 0 : NUkl[2][2] = facMS
1293 0 : * conj(coupSUSYPtr->getRsqqX(iGen3,idIn2A,k))
1294 0 : * conj(coupSUSYPtr->getRsqqX(iGen4,idIn1A,k))
1295 0 : * coupSUSYPtr->getRsqqX(iGen3,idIn2A,l)
1296 0 : * coupSUSYPtr->getRsqqX(iGen4,idIn1A,l);
1297 0 : NTUkl[1][1] = facMS
1298 0 : * real( conj(coupSUSYPtr->getLsqqX(iGen4,idIn2A,k))
1299 0 : * conj(coupSUSYPtr->getLsqqX(iGen3,idIn1A,k))
1300 0 : * coupSUSYPtr->getLsqqX(iGen3,idIn2A,l)
1301 0 : * coupSUSYPtr->getLsqqX(iGen4,idIn1A,l) );
1302 0 : NTUkl[1][2] = facTU
1303 0 : * real( conj(coupSUSYPtr->getRsqqX(iGen4,idIn2A,k))
1304 0 : * conj(coupSUSYPtr->getLsqqX(iGen3,idIn1A,k))
1305 0 : * coupSUSYPtr->getRsqqX(iGen3,idIn2A,l)
1306 0 : * coupSUSYPtr->getLsqqX(iGen4,idIn1A,l) );
1307 0 : NTUkl[2][1] = facTU
1308 0 : * real( conj(coupSUSYPtr->getLsqqX(iGen4,idIn2A,k))
1309 0 : * conj(coupSUSYPtr->getRsqqX(iGen3,idIn1A,k))
1310 0 : * coupSUSYPtr->getLsqqX(iGen3,idIn2A,l)
1311 0 : * coupSUSYPtr->getRsqqX(iGen4,idIn1A,l) );
1312 0 : NTUkl[2][2] = facMS
1313 0 : * real( conj(coupSUSYPtr->getRsqqX(iGen4,idIn2A,k))
1314 0 : * conj(coupSUSYPtr->getRsqqX(iGen3,idIn1A,k))
1315 0 : * coupSUSYPtr->getRsqqX(iGen3,idIn2A,l)
1316 0 : * coupSUSYPtr->getRsqqX(iGen4,idIn1A,l) );
1317 :
1318 : // Add to sums
1319 0 : sumNt += sigmaNeut / tNeut[k] / tNeut[l]
1320 0 : * real(NTkl[1][1] + NTkl[1][2] + NTkl[2][1] + NTkl[2][2]);
1321 0 : sumNu += sigmaNeut / uNeut[k] / uNeut[l]
1322 0 : * real(NUkl[1][1] + NUkl[1][2] + NUkl[2][1] + NUkl[2][2]);
1323 0 : sumInterference += 2.0 / 3.0 * sigmaNeut
1324 0 : * real(NTUkl[1][1] + NTUkl[1][2] + NTUkl[2][1] + NTUkl[2][2])
1325 0 : / tNeut[k] / uNeut[l];
1326 0 : }
1327 :
1328 : // Neutralino / Gluino interference
1329 :
1330 : // k-dependent factor for LL and RR contributions
1331 0 : double facMS = sH * particleDataPtr->m0(coupSUSYPtr->idNeut(k))
1332 0 : * particleDataPtr->m0(1000021);
1333 :
1334 : // Note: Nxkl defined as in [Boz07] with sigmaNeutGlu factored out
1335 : // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
1336 0 : complex NGA[3][3], NGB[3][3];
1337 0 : NGA[1][1] = facMS
1338 0 : * real( conj(coupSUSYPtr->getLsqqX(iGen4,idIn2A,k))
1339 0 : * conj(coupSUSYPtr->getLsqqX(iGen3,idIn1A,k))
1340 0 : * conj(coupSUSYPtr->getLsqqG(iGen3,idIn2A))
1341 0 : * conj(coupSUSYPtr->getLsqqG(iGen4,idIn1A)) );
1342 0 : NGA[1][2] = facTU
1343 0 : * real( conj(coupSUSYPtr->getRsqqX(iGen4,idIn2A,k))
1344 0 : * conj(coupSUSYPtr->getLsqqX(iGen3,idIn1A,k))
1345 0 : * conj(coupSUSYPtr->getLsqqG(iGen3,idIn2A))
1346 0 : * conj(coupSUSYPtr->getRsqqG(iGen4,idIn1A)) );
1347 0 : NGA[2][1] = facTU
1348 0 : * real( conj(coupSUSYPtr->getLsqqX(iGen4,idIn2A,k))
1349 0 : * conj(coupSUSYPtr->getRsqqX(iGen3,idIn1A,k))
1350 0 : * conj(coupSUSYPtr->getRsqqG(iGen3,idIn2A))
1351 0 : * conj(coupSUSYPtr->getLsqqG(iGen4,idIn1A)) );
1352 0 : NGA[2][2] = facMS
1353 0 : * real( conj(coupSUSYPtr->getRsqqX(iGen4,idIn2A,k))
1354 0 : * conj(coupSUSYPtr->getRsqqX(iGen3,idIn1A,k))
1355 0 : * conj(coupSUSYPtr->getRsqqG(iGen3,idIn2A))
1356 0 : * conj(coupSUSYPtr->getRsqqG(iGen4,idIn1A)) );
1357 0 : NGB[1][1] = facMS
1358 0 : * real( conj(coupSUSYPtr->getLsqqX(iGen3,idIn2A,k))
1359 0 : * conj(coupSUSYPtr->getLsqqX(iGen4,idIn1A,k))
1360 0 : * conj(coupSUSYPtr->getLsqqG(iGen4,idIn2A))
1361 0 : * conj(coupSUSYPtr->getLsqqG(iGen3,idIn1A)) );
1362 0 : NGB[1][2] = facMS
1363 0 : * real( conj(coupSUSYPtr->getRsqqX(iGen3,idIn2A,k))
1364 0 : * conj(coupSUSYPtr->getLsqqX(iGen4,idIn1A,k))
1365 0 : * conj(coupSUSYPtr->getRsqqG(iGen4,idIn2A))
1366 0 : * conj(coupSUSYPtr->getLsqqG(iGen3,idIn1A)) );
1367 0 : NGB[2][1] = facMS
1368 0 : * real( conj(coupSUSYPtr->getLsqqX(iGen3,idIn2A,k))
1369 0 : * conj(coupSUSYPtr->getRsqqX(iGen4,idIn1A,k))
1370 0 : * conj(coupSUSYPtr->getLsqqG(iGen4,idIn2A))
1371 0 : * conj(coupSUSYPtr->getRsqqG(iGen3,idIn1A)) );
1372 0 : NGB[2][2] = facMS
1373 0 : * real( conj(coupSUSYPtr->getRsqqX(iGen3,idIn2A,k))
1374 0 : * conj(coupSUSYPtr->getRsqqX(iGen4,idIn1A,k))
1375 0 : * conj(coupSUSYPtr->getRsqqG(iGen4,idIn2A))
1376 0 : * conj(coupSUSYPtr->getRsqqG(iGen3,idIn1A)) );
1377 :
1378 : // Add to sums
1379 0 : sumInterference += sigmaNeutGlu *
1380 0 : ( real(NGA[1][1] + NGA[1][2] + NGA[2][1] + NGA[2][2])
1381 0 : / tNeut[k] / uGlu
1382 0 : + real(NGB[1][1] + NGB[1][2] + NGB[2][1] + NGB[2][2])
1383 0 : / uNeut[k] / tGlu );
1384 0 : }
1385 0 : }
1386 : // t-channel + u-channel Gluinos + t/u interference
1387 :
1388 : // factor for LL and RR contributions
1389 0 : double facMS = sH * m2Glu;
1390 :
1391 : // Note: GT, GU defined as in [Boz07] with sigmaGlu factored out
1392 : // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
1393 0 : complex GT[3][3], GU[3][3], GTU[3][3];
1394 0 : GT[1][1] = facMS
1395 0 : * norm(coupSUSYPtr->getLsqqG(iGen4,idIn2A)
1396 0 : * coupSUSYPtr->getLsqqG(iGen3,idIn1A));
1397 0 : GT[1][2] = facTU
1398 0 : * norm(coupSUSYPtr->getRsqqG(iGen4,idIn2A)
1399 0 : * coupSUSYPtr->getLsqqG(iGen3,idIn1A));
1400 0 : GT[2][1] = facTU
1401 0 : * norm(coupSUSYPtr->getLsqqG(iGen4,idIn2A)
1402 0 : * coupSUSYPtr->getRsqqG(iGen3,idIn1A));
1403 0 : GT[2][2] = facMS
1404 0 : * norm(coupSUSYPtr->getRsqqG(iGen4,idIn2A)
1405 0 : * coupSUSYPtr->getRsqqG(iGen3,idIn1A));
1406 0 : GU[1][1] = facMS
1407 0 : * norm(coupSUSYPtr->getLsqqG(iGen3,idIn2A)
1408 0 : * coupSUSYPtr->getLsqqG(iGen4,idIn1A));
1409 0 : GU[1][2] = facTU
1410 0 : * norm(coupSUSYPtr->getLsqqG(iGen3,idIn2A)
1411 0 : * coupSUSYPtr->getRsqqG(iGen4,idIn1A));
1412 0 : GU[2][1] = facTU
1413 0 : * norm(coupSUSYPtr->getRsqqG(iGen3,idIn2A)
1414 0 : * coupSUSYPtr->getLsqqG(iGen4,idIn1A));
1415 0 : GU[2][2] = facMS
1416 0 : * norm(coupSUSYPtr->getRsqqG(iGen3,idIn2A)
1417 0 : * coupSUSYPtr->getRsqqG(iGen4,idIn1A));
1418 :
1419 0 : GTU[1][1] = facMS
1420 0 : * real(coupSUSYPtr->getLsqqG(iGen3,idIn1A)
1421 0 : * coupSUSYPtr->getLsqqG(iGen4,idIn2A)
1422 0 : * conj(coupSUSYPtr->getLsqqG(iGen3,idIn2A))
1423 0 : * conj(coupSUSYPtr->getLsqqG(iGen4,idIn1A)) );
1424 :
1425 0 : GTU[1][2] = facTU
1426 0 : * real(coupSUSYPtr->getLsqqG(iGen3,idIn1A)
1427 0 : * coupSUSYPtr->getRsqqG(iGen4,idIn2A)
1428 0 : * conj(coupSUSYPtr->getRsqqG(iGen3,idIn2A))
1429 0 : * conj(coupSUSYPtr->getLsqqG(iGen4,idIn1A)) );
1430 :
1431 0 : GTU[2][1] = facTU
1432 0 : * real(coupSUSYPtr->getRsqqG(iGen3,idIn1A)
1433 0 : * coupSUSYPtr->getLsqqG(iGen4,idIn2A)
1434 0 : * conj(coupSUSYPtr->getLsqqG(iGen3,idIn2A))
1435 0 : * conj(coupSUSYPtr->getRsqqG(iGen4,idIn1A)) );
1436 :
1437 0 : GTU[2][2] = facMS
1438 0 : * real(coupSUSYPtr->getRsqqG(iGen3,idIn1A)
1439 0 : * coupSUSYPtr->getRsqqG(iGen4,idIn2A)
1440 0 : * conj(coupSUSYPtr->getRsqqG(iGen3,idIn2A))
1441 0 : * conj(coupSUSYPtr->getRsqqG(iGen4,idIn1A)) );
1442 :
1443 : // Add to sums
1444 0 : sumGt += sigmaGlu * real(GT[1][1] + GT[1][2] + GT[2][1] + GT[2][2])
1445 0 : / pow2(tGlu) ;
1446 0 : sumGu += sigmaGlu * real(GU[1][1] + GU[1][2] + GU[2][1] + GU[2][2])
1447 0 : / pow2(uGlu) ;
1448 0 : sumInterference += - 2.0 / 3.0 * sigmaGlu
1449 0 : * real(GTU[1][1] + GTU[1][2] + GTU[2][1] + GTU[2][2])
1450 0 : / tGlu / uGlu;
1451 :
1452 0 : }
1453 :
1454 : // Cross section
1455 0 : double sigma = sumNt + sumNu + sumCt + sumCu + sumGt + sumGu
1456 0 : + sumInterference;
1457 :
1458 : // Identical particles?
1459 0 : if (id3Sav == id4Sav) sigma /= 2.0;
1460 :
1461 : // Return answer.
1462 : return sigma;
1463 :
1464 0 : }
1465 :
1466 : //--------------------------------------------------------------------------
1467 :
1468 : // Select identity, colour and anticolour.
1469 :
1470 : void Sigma2qq2squarksquark::setIdColAcol() {
1471 :
1472 : // Set flavours.
1473 0 : if (id1 > 0 && id2 > 0) {
1474 0 : setId( id1, id2, id3Sav, id4Sav);
1475 0 : } else {
1476 : // 1,2 -> -3,-4
1477 0 : setId( id1, id2,-id3Sav,-id4Sav);
1478 : }
1479 :
1480 : // Coded sigma is for ud -> ~q~q'. Swap t and u for du -> ~q~q'.
1481 0 : swapTU = (isUD && abs(id1) % 2 == 0);
1482 :
1483 : // Select colour flow topology
1484 : // Recompute contributions to this particular in- out- flavour combination
1485 0 : sigmaHat();
1486 : // A: t-channel neutralino, t-channel chargino, or u-channel gluino
1487 0 : double sumA = sumNt + sumCt + sumGu;
1488 0 : double sumAB = sumNt + sumNu + sumCt + sumCu + sumGt + sumGu;
1489 0 : if (swapTU) sumA = sumAB - sumA;
1490 0 : setColAcol( 1, 0, 2, 0, 1, 0, 2, 0);
1491 : // B: t-channel gluino or u-channel neutralino
1492 0 : if (rndmPtr->flat()*sumAB > sumA) setColAcol( 1, 0, 2, 0, 2, 0, 1, 0);
1493 :
1494 : // Switch to anti-colors if antiquarks
1495 0 : if (id1 < 0 || id2 < 0) swapColAcol();
1496 :
1497 0 : }
1498 :
1499 : //==========================================================================
1500 :
1501 : // Sigma2qqbar2squarkantisquark
1502 : // Cross section for qqbar-initiated squark-antisquark production
1503 :
1504 : //--------------------------------------------------------------------------
1505 :
1506 : // Initialize process.
1507 :
1508 : void Sigma2qqbar2squarkantisquark::initProc() {
1509 :
1510 : //Typecast to the correct couplings
1511 0 : coupSUSYPtr = (CoupSUSY*) couplingsPtr;
1512 :
1513 : // Is this a ~u_i ~d*_j, ~d_i ~u*_j final state or ~d_i ~d*_j, ~u_i ~u*_j
1514 0 : if (abs(id3Sav) % 2 == abs(id4Sav) % 2) isUD = false;
1515 0 : else isUD = true;
1516 :
1517 : // Extract isospin and mass-ordering indices
1518 0 : if(isUD && abs(id3Sav)%2 == 1){
1519 0 : iGen3 = 3*(abs(id4Sav)/2000000) + (abs(id3Sav)%10+1)/2;
1520 0 : iGen4 = 3*(abs(id3Sav)/2000000) + (abs(id4Sav)%10+1)/2;
1521 0 : }
1522 : else {
1523 0 : iGen3 = 3*(abs(id3Sav)/2000000) + (abs(id3Sav)%10+1)/2;
1524 0 : iGen4 = 3*(abs(id4Sav)/2000000) + (abs(id4Sav)%10+1)/2;
1525 : }
1526 :
1527 : // Derive name
1528 0 : nameSave = "q qbar' -> "+particleDataPtr->name(abs(id3Sav))+" "+
1529 0 : particleDataPtr->name(-abs(id4Sav));
1530 0 : if (isUD && abs(id3Sav) != abs(id4Sav)) nameSave +=" + c.c.";
1531 :
1532 : // Count 5 neutralinos in NMSSM
1533 0 : nNeut = (coupSUSYPtr->isNMSSM ? 5 : 4);
1534 :
1535 : // Store mass squares of all possible internal propagator lines
1536 0 : m2Glu = pow2(particleDataPtr->m0(1000021));
1537 0 : m2Neut.resize(nNeut+1);
1538 0 : for (int iNeut=1;iNeut<=nNeut;iNeut++)
1539 0 : m2Neut[iNeut] = pow2(particleDataPtr->m0(coupSUSYPtr->idNeut(iNeut)));
1540 :
1541 : // Set sizes of some arrays to be used below
1542 0 : tNeut.resize(nNeut+1);
1543 0 : uNeut.resize(nNeut+1);
1544 :
1545 : // Shorthand for Weak mixing
1546 0 : xW = coupSUSYPtr->sin2W;
1547 :
1548 : // Secondary open width fraction.
1549 0 : openFracPair = particleDataPtr->resOpenFrac(id3Sav, id4Sav);
1550 :
1551 : // Select interference terms
1552 0 : onlyQCD = settingsPtr->flag("SUSY:qqbar2squarkantisquark:onlyQCD");
1553 0 : }
1554 :
1555 : //--------------------------------------------------------------------------
1556 :
1557 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
1558 :
1559 : void Sigma2qqbar2squarkantisquark::sigmaKin() {
1560 :
1561 : // Z/W propagator
1562 0 : if (! isUD) {
1563 0 : double sV= sH - pow2(coupSUSYPtr->mZpole);
1564 0 : double d = pow2(sV) + pow2(coupSUSYPtr->mZpole * coupSUSYPtr->wZpole);
1565 0 : propZW = complex( sV / d, coupSUSYPtr->mZpole * coupSUSYPtr->wZpole / d);
1566 0 : } else {
1567 0 : double sV= sH - pow2(coupSUSYPtr->mWpole);
1568 0 : double d = pow2(sV) + pow2(coupSUSYPtr->mWpole * coupSUSYPtr->wWpole);
1569 0 : propZW = complex( sV / d, coupSUSYPtr->mWpole * coupSUSYPtr->wWpole / d);
1570 0 : }
1571 :
1572 : // Flavor-independent pre-factors
1573 0 : double comFacHat = M_PI/sH2 * openFracPair;
1574 :
1575 0 : sigmaEW = comFacHat * pow2(alpEM);
1576 0 : sigmaGlu = comFacHat * 2.0 * pow2(alpS) / 9.0;
1577 0 : sigmaEWG = comFacHat * 8.0 * alpEM * alpS / 9.0;
1578 :
1579 0 : }
1580 :
1581 : //--------------------------------------------------------------------------
1582 :
1583 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1584 :
1585 : double Sigma2qqbar2squarkantisquark::sigmaHat() {
1586 :
1587 : // In-pair must be opposite-sign
1588 0 : if (id1 * id2 > 0) return 0.0;
1589 :
1590 : // Check correct charge sum
1591 0 : if (isUD && abs(id1) %2 == abs(id2) % 2) return 0.0;
1592 0 : if (!isUD && abs(id1) % 2 != abs(id2) % 2) return 0.0;
1593 :
1594 : // Check if using QCD diagrams only
1595 :
1596 :
1597 : // Coded UD sigma is for udbar -> ~u~d'*. Swap t<->u for dbaru -> ~u~d'*.
1598 0 : swapTU = (isUD && abs(id1) % 2 != 0);
1599 :
1600 : // Coded QQ sigma is for qqbar -> ~q~q*. Swap t<->u for qbarq -> ~q~q*.
1601 0 : if (!isUD && id1 < 0) swapTU = true;
1602 :
1603 : // Generation indices of incoming particles
1604 0 : int idIn1A = (swapTU) ? abs(id2) : abs(id1);
1605 0 : int idIn2A = (swapTU) ? abs(id1) : abs(id2);
1606 0 : int iGen1 = (idIn1A+1)/2;
1607 0 : int iGen2 = (idIn2A+1)/2;
1608 :
1609 : // Auxiliary factors for use below
1610 0 : tGlu = tH - m2Glu;
1611 0 : uGlu = uH - m2Glu;
1612 0 : for (int i=1; i<= nNeut; i++) {
1613 0 : tNeut[i] = tH - m2Neut[i];
1614 0 : uNeut[i] = uH - m2Neut[i];
1615 : }
1616 :
1617 : // Initial values for pieces used for color-flow selection below
1618 0 : sumColS = 0.0;
1619 0 : sumColT = 0.0;
1620 0 : sumInterference = 0.0;
1621 :
1622 : // Common factor for LR and RL contributions
1623 0 : double facTU = uH*tH-s3*s4;
1624 :
1625 : // Case A) Opposite-isospin: udbar -> ~u~d*
1626 0 : if ( isUD ) {
1627 :
1628 : // s-channel W contribution (only contributes to LL helicities)
1629 0 : if ( !onlyQCD ) {
1630 0 : sumColS += sigmaEW / 16.0 / pow2(xW) / pow2(1.0-xW)
1631 0 : * norm(conj(coupSUSYPtr->LudW[iGen1][iGen2])
1632 0 : * coupSUSYPtr->LsusdW[iGen3][iGen4]) * facTU
1633 0 : * norm(propZW);
1634 0 : }
1635 :
1636 : // t-channel gluino contributions
1637 : double GT[3][3];
1638 0 : double facLR = m2Glu * sH;
1639 : // LL, LR, RL, RR
1640 0 : GT[1][1] = facTU * norm(conj(coupSUSYPtr->LsddG[iGen4][iGen2])
1641 0 : *coupSUSYPtr->LsuuG[iGen3][iGen1]);
1642 0 : GT[1][2] = facLR * norm(conj(coupSUSYPtr->RsddG[iGen4][iGen2])
1643 0 : *coupSUSYPtr->LsuuG[iGen3][iGen1]);
1644 0 : GT[2][1] = facLR * norm(conj(coupSUSYPtr->LsddG[iGen4][iGen2])
1645 0 : *coupSUSYPtr->RsuuG[iGen3][iGen1]);
1646 0 : GT[2][2] = facTU * norm(conj(coupSUSYPtr->RsddG[iGen4][iGen2])
1647 0 : *coupSUSYPtr->RsuuG[iGen3][iGen1]);
1648 : // leading color flow for t-channel gluino is annihilation-like
1649 0 : sumColS += sigmaGlu / pow2(tGlu)
1650 0 : * (GT[1][1] + GT[1][2] + GT[2][1] + GT[2][2]);
1651 :
1652 : // W-Gluino interference (only contributes to LL helicities)
1653 0 : if ( !onlyQCD ) {
1654 0 : sumColS += sigmaEWG / 4.0 / xW / (1-xW)
1655 0 : * real(conj(coupSUSYPtr->LsuuG[iGen3][iGen1])
1656 0 : * coupSUSYPtr->LsddG[iGen4][iGen2]
1657 0 : * conj(coupSUSYPtr->LudW[iGen1][iGen2])
1658 0 : * coupSUSYPtr->LsusdW[iGen3][iGen4]) * facTU
1659 0 : / tGlu * sqrt(norm(propZW));
1660 0 : }
1661 :
1662 : // t-channel neutralinos
1663 : // NOT YET IMPLEMENTED !
1664 :
1665 0 : }
1666 :
1667 : // Case B) Same-isospin: qqbar -> ~d~d* , ~u~u*
1668 : else {
1669 :
1670 0 : double eQ = (idIn1A % 2 == 0) ? 2./3. : 1./3. ;
1671 0 : double eSq = (abs(id3Sav) % 2 == 0) ? 2./3. : 1./3. ;
1672 :
1673 : // s-channel gluon (strictly flavor-diagonal)
1674 0 : if (abs(id3Sav) == abs(id4Sav) && abs(id1) == abs(id2)) {
1675 : // Factor 2 since contributes to both ha != hb helicities
1676 0 : sumColT += 2. * sigmaGlu * facTU / pow2(sH);
1677 0 : }
1678 :
1679 : // t-channel gluino (only for in-isospin = out-isospin).
1680 0 : if (eQ == eSq) {
1681 : // Sum over helicities.
1682 : double GT[3][3];
1683 0 : double facLR = sH * m2Glu;
1684 0 : GT[1][1] = facTU * norm(coupSUSYPtr->getLsqqG(iGen3,idIn1A)
1685 0 : * conj(coupSUSYPtr->getLsqqG(iGen4,idIn2A)));
1686 0 : GT[1][2] = facLR * norm(coupSUSYPtr->getLsqqG(iGen3,idIn1A)
1687 0 : * conj(coupSUSYPtr->getRsqqG(iGen4,idIn2A)));
1688 0 : GT[2][1] = facLR * norm(coupSUSYPtr->getRsqqG(iGen3,idIn1A)
1689 0 : * conj(coupSUSYPtr->getLsqqG(iGen4,idIn2A)));
1690 0 : GT[2][2] = facTU * norm(coupSUSYPtr->getRsqqG(iGen3,idIn1A)
1691 0 : * conj(coupSUSYPtr->getRsqqG(iGen4,idIn2A)));
1692 : // Add contribution to color topology: S
1693 0 : sumColS += sigmaGlu / pow2(tGlu)
1694 0 : * ( GT[1][1] + GT[1][2] + GT[2][1] + GT[2][2]);
1695 :
1696 : // gluon-gluino interference (strictly flavor-diagonal)
1697 0 : if (abs(id3Sav) == abs(id4Sav) && abs(id1) == abs(id2)) {
1698 : double GG11, GG22;
1699 0 : GG11 = - facTU * 2./3.
1700 0 : * real( conj(coupSUSYPtr->getLsqqG(iGen3,idIn1A))
1701 0 : * coupSUSYPtr->getLsqqG(iGen4,idIn2A));
1702 : GG22 = - facTU * 2./3.
1703 0 : * real( conj(coupSUSYPtr->getRsqqG(iGen3,idIn1A))
1704 0 : * coupSUSYPtr->getRsqqG(iGen4,idIn2A));
1705 : // Sum over two contributing helicities
1706 0 : sumInterference += sigmaGlu / sH / tGlu
1707 0 : * ( GG11 + GG22 );
1708 0 : }
1709 :
1710 0 : }
1711 :
1712 : // Skip the rest if only including QCD diagrams
1713 0 : if (onlyQCD) return sumColT+sumColS+sumInterference;
1714 :
1715 : // s-channel photon (strictly flavor-diagonal) and Z/gamma interference
1716 0 : if (abs(id3Sav) == abs(id4Sav) && abs(id1) == abs(id2)) {
1717 :
1718 : // gamma
1719 : // Factor 2 since contributes to both ha != hb helicities
1720 0 : sumColS += 2. * pow2(eQ) * pow2(eSq) * sigmaEW * facTU / pow2(sH);
1721 :
1722 : // Z/gamma interference
1723 0 : double CsqZ = real(coupSUSYPtr->LsusuZ[iGen3][iGen4]
1724 0 : + coupSUSYPtr->RsusuZ[iGen3][iGen4]);
1725 0 : if (abs(id3Sav)%2 != 0) CsqZ = real(coupSUSYPtr->LsdsdZ[iGen3][iGen4]
1726 0 : + coupSUSYPtr->RsdsdZ[iGen3][iGen4]);
1727 0 : sumColS += eQ * eSq * sigmaEW * facTU / 2.0 / xW / (1.-xW)
1728 0 : * sqrt(norm(propZW)) / sH * CsqZ
1729 0 : * (coupSUSYPtr->LqqZ[idIn1A] + coupSUSYPtr->LqqZ[idIn2A]);
1730 :
1731 : // Gluino/gamma interference (only for same-isospin)
1732 0 : if (eQ == eSq) {
1733 0 : double CsqG11 = real(conj(coupSUSYPtr->LsuuG[iGen3][iGen1])
1734 0 : *coupSUSYPtr->LsuuG[iGen4][iGen2]);
1735 0 : double CsqG22 = real(conj(coupSUSYPtr->RsuuG[iGen3][iGen1])
1736 0 : *coupSUSYPtr->RsuuG[iGen4][iGen2]);
1737 0 : if (id3Sav%2 != 0) {
1738 0 : CsqG11 = real(conj(coupSUSYPtr->LsddG[iGen3][iGen1])
1739 0 : *coupSUSYPtr->LsddG[iGen4][iGen2]);
1740 0 : CsqG22 = real(conj(coupSUSYPtr->RsddG[iGen3][iGen1])
1741 0 : *coupSUSYPtr->RsddG[iGen4][iGen2]);
1742 0 : }
1743 0 : sumColS += eQ * eSq * sigmaEWG * facTU
1744 0 : * (CsqG11 + CsqG22) / sH / tGlu;
1745 0 : }
1746 0 : }
1747 :
1748 : // s-channel Z (only for q flavor = qbar flavor)
1749 0 : if (abs(id1) == abs(id2)) {
1750 0 : double CsqZ = norm(coupSUSYPtr->LsusuZ[iGen3][iGen4]
1751 0 : + coupSUSYPtr->RsusuZ[iGen3][iGen4]);
1752 0 : if (abs(id3Sav)%2 != 0) CsqZ = norm(coupSUSYPtr->LsdsdZ[iGen3][iGen4]
1753 0 : + coupSUSYPtr->RsdsdZ[iGen3][iGen4]);
1754 0 : sumColS += sigmaEW * facTU / 16.0 / pow2(xW) / pow2(1.0-xW)
1755 0 : * norm(propZW) * CsqZ * ( pow2(coupSUSYPtr->LqqZ[idIn1A])
1756 0 : + pow2(coupSUSYPtr->RqqZ[idIn1A]) );
1757 :
1758 : // Z/gluino interference (only for in-isospin = out-isospin)
1759 0 : if (eQ == eSq) {
1760 0 : double GZ11 = real(conj(coupSUSYPtr->getLsqqG(iGen3,idIn1A))
1761 0 : *coupSUSYPtr->getLsqqG(iGen4,idIn2A)
1762 0 : *(coupSUSYPtr->getLsqsqZ(id3Sav,id4Sav)
1763 0 : +coupSUSYPtr->getRsqsqZ(id3Sav,id4Sav)))
1764 0 : *coupSUSYPtr->LqqZ[idIn1A];
1765 0 : double GZ22 = real(conj(coupSUSYPtr->getRsqqG(iGen3,idIn1A))
1766 0 : *coupSUSYPtr->getRsqqG(iGen4,idIn2A)
1767 0 : *(coupSUSYPtr->getLsqsqZ(id3Sav,id4Sav)
1768 0 : +coupSUSYPtr->getRsqsqZ(id3Sav,id4Sav)))
1769 0 : *coupSUSYPtr->RqqZ[idIn1A];
1770 0 : sumColS += sigmaEWG * facTU / 4.0 / xW / (1.-xW)
1771 0 : * ( GZ11 + GZ22 ) * sqrt(norm(propZW)) / tGlu;
1772 0 : }
1773 0 : }
1774 :
1775 : // t-channel neutralinos
1776 : // NOT YET IMPLEMENTED !
1777 :
1778 0 : }
1779 :
1780 : // Cross section
1781 0 : double sigma = sumColS + sumColT + sumInterference;
1782 :
1783 : // Return answer.
1784 : return sigma;
1785 :
1786 0 : }
1787 :
1788 : //--------------------------------------------------------------------------
1789 :
1790 : // Select identity, colour and anticolour.
1791 :
1792 : void Sigma2qqbar2squarkantisquark::setIdColAcol() {
1793 :
1794 : // Check if charge conjugate final state?
1795 0 : isCC = false;
1796 0 : if (isUD && ( (id1-1)%2 < 0 || (id2-1)%2 < 0 )) isCC = true;
1797 :
1798 : //check if charge conjugate
1799 0 : id3 = (isCC) ? -id3Sav : id3Sav;
1800 0 : id4 = (isCC) ? -id4Sav : id4Sav;
1801 :
1802 : // Set flavours.
1803 0 : setId( id1, id2, id3, id4);
1804 :
1805 : // Coded UD sigma is for udbar -> ~u~d'*. Swap t<->u for dbaru -> ~u~d'*.
1806 : // Coded QQ sigma is for qqbar -> ~q~q*. Swap t<->u for qbarq -> ~q~q*.
1807 0 : if (isUD) {
1808 0 : swapTU = (abs(id1) % 2 != 0);
1809 0 : } else {
1810 0 : swapTU = (id1 < 0);
1811 : }
1812 :
1813 : // Select colour flow topology
1814 : // Recompute individual contributions to this in-out flavour combination
1815 0 : sigmaHat();
1816 0 : double R = rndmPtr->flat();
1817 0 : double fracS = sumColS / (sumColS + sumColT) ;
1818 : // S: color flow as in S-channel singlet
1819 0 : if (R < fracS) {
1820 0 : setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
1821 0 : if (swapTU) setColAcol( 0, 1, 1, 0, 2, 0, 0, 2);
1822 : }
1823 : // T: color flow as in T-channel singlet
1824 : else {
1825 0 : setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
1826 0 : if (swapTU) setColAcol( 0, 1, 2, 0, 2, 0, 0, 1);
1827 : }
1828 :
1829 0 : if (isCC) swapColAcol();
1830 :
1831 0 : }
1832 :
1833 : //==========================================================================
1834 :
1835 : // Sigma2gg2squarkantisquark
1836 : // Cross section for gg-initiated squark-antisquark production
1837 :
1838 : //--------------------------------------------------------------------------
1839 :
1840 : // Initialize process.
1841 :
1842 : void Sigma2gg2squarkantisquark::initProc() {
1843 :
1844 : //Typecast to the correct couplings
1845 0 : coupSUSYPtr = (CoupSUSY*) couplingsPtr;
1846 :
1847 : // Process Name
1848 0 : nameSave = "g g -> "+particleDataPtr->name(abs(id3Sav))+" "
1849 0 : +particleDataPtr->name(-abs(id4Sav));
1850 :
1851 : // Squark pole mass
1852 0 : m2Sq = pow2(particleDataPtr->m0(id3Sav));
1853 :
1854 : // Secondary open width fraction.
1855 0 : openFracPair = particleDataPtr->resOpenFrac(id3Sav, id4Sav);
1856 :
1857 0 : }
1858 :
1859 : //--------------------------------------------------------------------------
1860 :
1861 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
1862 :
1863 : void Sigma2gg2squarkantisquark::sigmaKin() {
1864 :
1865 : // Modified Mandelstam variables for massive kinematics with m3 = m4.
1866 : // tHSq = tHat - m_squark^2; uHSq = uHat - m_squark^2.
1867 : // double s34Avg = 0.5 * (s3 + s4) - 0.25 * pow2(s3 - s4) / sH;
1868 0 : double tHSq = -0.5 * (sH - tH + uH);
1869 0 : double uHSq = -0.5 * (sH + tH - uH);
1870 : // ! (NEED TO CHECK THAT THESE APPLIED CORRECTLY BELOW) !
1871 : // ! (PRELIMINARY CROSS-CHECKS WITH PYTHIA 6 COME OUT OK) !
1872 :
1873 : // Helicity-independent prefactor
1874 0 : double comFacHat = M_PI/sH2 * pow2(alpS) / 128.0
1875 0 : * ( 24.0 * (1.0 - 2*tHSq*uHSq/sH2) - 8.0/3.0 );
1876 :
1877 : // Helicity-dependent factors
1878 0 : sigma = 0.0;
1879 0 : for (int ha=-1;ha<=1;ha += 2) {
1880 0 : for (int hb=-1;hb<=1;hb += 2) {
1881 : // Divide by 4 for helicity average
1882 0 : sigma += comFacHat / 4.0
1883 0 : * ( (1.0-ha*hb)
1884 0 : - 2.0 * sH*m2Sq/tHSq/uHSq
1885 0 : * ( 1.0 - ha*hb - sH*m2Sq/tHSq/uHSq));
1886 : }
1887 : }
1888 :
1889 0 : }
1890 :
1891 : //--------------------------------------------------------------------------
1892 :
1893 : // Select identity, colour and anticolour.
1894 :
1895 : void Sigma2gg2squarkantisquark::setIdColAcol() {
1896 :
1897 : // Set flavours.
1898 0 : setId( id1, id2, id3Sav, id4Sav);
1899 :
1900 : // Set color flow (random for now)
1901 0 : double R = rndmPtr->flat();
1902 0 : if (R < 0.5) setColAcol( 1, 2, 2, 3, 1, 0, 0, 3);
1903 0 : else setColAcol( 1, 2, 3, 1, 3, 0, 0, 2);
1904 :
1905 0 : }
1906 :
1907 : //==========================================================================
1908 :
1909 : // Sigma2qg2squarkgluino
1910 : // Cross section for squark-gluino production
1911 :
1912 : //--------------------------------------------------------------------------
1913 :
1914 : // Initialize process.
1915 :
1916 : void Sigma2qg2squarkgluino::initProc() {
1917 :
1918 : //Typecast to the correct couplings
1919 0 : coupSUSYPtr = (CoupSUSY*) couplingsPtr;
1920 :
1921 : // Derive name
1922 0 : nameSave = "q g -> "+particleDataPtr->name(abs(id3Sav))+" gluino + c.c.";
1923 :
1924 : // Final-state mass squares
1925 0 : m2Glu = pow2(particleDataPtr->m0(1000021));
1926 0 : m2Sq = pow2(particleDataPtr->m0(id3Sav));
1927 :
1928 : // Secondary open width fraction.
1929 0 : openFracPair = particleDataPtr->resOpenFrac(id3Sav, 1000021);
1930 :
1931 0 : }
1932 :
1933 : //--------------------------------------------------------------------------
1934 :
1935 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
1936 :
1937 : void Sigma2qg2squarkgluino::sigmaKin() {
1938 :
1939 : // Common pre-factor
1940 0 : comFacHat = (M_PI / sH2) * pow2(alpS) * 0.5 * openFracPair;
1941 :
1942 : // Invariants (still with Pythia 6 sign convention)
1943 0 : double tGlu = m2Glu-tH;
1944 0 : double uGlu = m2Glu-uH;
1945 0 : double tSq = m2Sq-tH;
1946 0 : double uSq = m2Sq-uH;
1947 :
1948 : // Color flow A: quark color annihilates with anticolor of g
1949 0 : sigmaA = 0.5*4./9.* tGlu/sH + (tGlu*sH+2.*m2Glu*tSq)/pow2(tGlu) -
1950 0 : ( (sH-m2Sq+m2Glu)*(-tSq)-sH*m2Glu )/sH/(-tGlu)
1951 0 : + 0.5*1./2.*( tSq*(tH+2.*uH+m2Glu)-tGlu*(sH-2.*tSq)
1952 0 : + (-uGlu)*(tH+m2Glu+2.*m2Sq) )/2./tGlu/uSq;
1953 : // Color flow B: quark and gluon colors iterchanged
1954 0 : sigmaB = 4./9.*(-uGlu)*(uH+m2Sq)/pow2(uSq)
1955 0 : + 1./18.* (sH*(uH+m2Glu) + 2.*(m2Sq-m2Glu)*uGlu)/sH/(-uSq)
1956 0 : + 0.5*4./9.*tGlu/sH
1957 0 : + 0.5*1./2.*(tSq*(tH+2.*uH+m2Glu)-tGlu*(sH-2.*tSq)
1958 0 : + (-uGlu)*(tH+m2Glu+2.*m2Sq))/2./tGlu/uSq;
1959 :
1960 0 : }
1961 :
1962 : //--------------------------------------------------------------------------
1963 :
1964 : double Sigma2qg2squarkgluino::sigmaHat() {
1965 :
1966 : // Check whether right incoming flavor
1967 0 : int idQA = (id1 == 21) ? id2 : id1;
1968 0 : int idSq = (abs(id3) == 10000021) ? id4 : id3;
1969 :
1970 : // Check for charge conservation
1971 0 : if(idQA%2 != idSq%2) return 0.0;
1972 :
1973 0 : int idQ = (abs(idQA)+1)/2;
1974 0 : idSq = 3 * (id3Sav / 2000000) + (id3Sav % 10 + 1)/2;
1975 :
1976 : double mixingFac;
1977 0 : if(abs(idQA) % 2 == 1)
1978 0 : mixingFac = norm(coupSUSYPtr->LsddG[idSq][idQ])
1979 0 : + norm(coupSUSYPtr->RsddG[idSq][idQ]);
1980 : else
1981 0 : mixingFac = norm(coupSUSYPtr->LsuuG[idSq][idQ])
1982 0 : + norm(coupSUSYPtr->RsuuG[idSq][idQ]);
1983 :
1984 0 : return mixingFac * comFacHat * (sigmaA + sigmaB);
1985 0 : }
1986 :
1987 : //--------------------------------------------------------------------------
1988 :
1989 : // Select identity, colour and anticolour.
1990 :
1991 : void Sigma2qg2squarkgluino::setIdColAcol() {
1992 :
1993 : // Check if charge conjugate final state?
1994 0 : int idQ = (id1 == 21) ? id2 : id1;
1995 0 : id3 = (idQ > 0) ? id3Sav : -id3Sav;
1996 0 : id4 = 1000021;
1997 :
1998 : // Set flavors
1999 0 : setId( id1, id2, id3, id4);
2000 :
2001 : // Select color flow A or B (see above)
2002 : // Recompute individual contributions to this in-out flavour combination
2003 0 : sigmaHat();
2004 0 : double R = rndmPtr->flat()*(sigmaA+sigmaB);
2005 0 : if (idQ == id1) {
2006 0 : setColAcol(1,0,2,1,3,0,2,3);
2007 0 : if (R > sigmaA) setColAcol(1,0,2,3,2,0,1,3);
2008 : } else {
2009 0 : setColAcol(2,1,1,0,3,0,2,3);
2010 0 : if (R > sigmaB) setColAcol(2,3,1,0,2,0,1,3);
2011 : }
2012 0 : if (idQ < 0) swapColAcol();
2013 :
2014 : // Use reflected kinematics if gq initial state
2015 0 : if (id1 == 21) swapTU = true;
2016 :
2017 0 : }
2018 :
2019 : //==========================================================================
2020 :
2021 : // Sigma2gg2gluinogluino
2022 : // Cross section for gluino pair production from gg initial states
2023 : // (validated against Pythia 6)
2024 :
2025 : //--------------------------------------------------------------------------
2026 :
2027 : // Initialize process.
2028 :
2029 : void Sigma2gg2gluinogluino::initProc() {
2030 :
2031 : //Typecast to the correct couplings
2032 0 : coupSUSYPtr = (CoupSUSY*) couplingsPtr;
2033 :
2034 : // Secondary open width fraction.
2035 0 : openFracPair = particleDataPtr->resOpenFrac(1000021, 1000021);
2036 :
2037 0 : }
2038 :
2039 : //--------------------------------------------------------------------------
2040 :
2041 : // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
2042 :
2043 : void Sigma2gg2gluinogluino::sigmaKin() {
2044 :
2045 : // Modified Mandelstam variables for massive kinematics with m3 = m4.
2046 : // tHG = tHat - m_gluino^2; uHG = uHat - m_gluino^2.
2047 0 : double s34Avg = 0.5 * (s3 + s4) - 0.25 * pow2(s3 - s4) / sH;
2048 0 : double tHG = -0.5 * (sH - tH + uH);
2049 0 : double uHG = -0.5 * (sH + tH - uH);
2050 0 : double tHG2 = tHG * tHG;
2051 0 : double uHG2 = uHG * uHG;
2052 :
2053 : // Calculate kinematics dependence.
2054 0 : sigTS = (tHG * uHG - 2. * s34Avg * (tHG + 2. * s34Avg)) / tHG2
2055 0 : + (tHG * uHG + s34Avg * (uHG - tHG)) / (sH * tHG);
2056 0 : sigUS = (tHG * uHG - 2. * s34Avg * (uHG + 2. * s34Avg)) / uHG2
2057 0 : + (tHG * uHG + s34Avg * (tHG - uHG)) / (sH * uHG);
2058 0 : sigTU = 2. * tHG * uHG / sH2 + s34Avg * (sH - 4. * s34Avg)
2059 0 : / (tHG * uHG);
2060 0 : sigSum = sigTS + sigUS + sigTU;
2061 :
2062 : // Answer contains factor 1/2 from identical gluinos.
2063 0 : sigma = (M_PI / sH2) * pow2(alpS) * (9./4.) * 0.5 * sigSum
2064 0 : * openFracPair;
2065 :
2066 0 : }
2067 :
2068 : //--------------------------------------------------------------------------
2069 :
2070 : // Select identity, colour and anticolour.
2071 :
2072 : void Sigma2gg2gluinogluino::setIdColAcol() {
2073 :
2074 : // Flavours are trivial.
2075 0 : setId( id1, id2, 1000021, 1000021);
2076 :
2077 : // Three colour flow topologies, each with two orientations.
2078 0 : double sigRand = sigSum * rndmPtr->flat();
2079 0 : if (sigRand < sigTS) setColAcol( 1, 2, 2, 3, 1, 4, 4, 3);
2080 0 : else if (sigRand < sigTS + sigUS)
2081 0 : setColAcol( 1, 2, 3, 1, 3, 4, 4, 2);
2082 0 : else setColAcol( 1, 2, 3, 4, 1, 4, 3, 2);
2083 0 : if (rndmPtr->flat() > 0.5) swapColAcol();
2084 :
2085 0 : }
2086 :
2087 : //==========================================================================
2088 :
2089 : // Sigma2qqbar2gluinogluino
2090 : // Cross section for gluino pair production from qqbar initial states
2091 : // (validated against Pythia 6 for SLHA1 case)
2092 :
2093 : //--------------------------------------------------------------------------
2094 :
2095 : // Initialize process.
2096 :
2097 : void Sigma2qqbar2gluinogluino::initProc() {
2098 :
2099 : //Typecast to the correct couplings
2100 0 : coupSUSYPtr = (CoupSUSY*) couplingsPtr;
2101 :
2102 : // Secondary open width fraction.
2103 0 : openFracPair = particleDataPtr->resOpenFrac(1000021, 1000021);
2104 :
2105 0 : }
2106 :
2107 : //--------------------------------------------------------------------------
2108 :
2109 : // Begin evaluate d(sigmaHat)/d(tHat); flavour-independent part.
2110 :
2111 : void Sigma2qqbar2gluinogluino::sigmaKin() {
2112 :
2113 : // Modified Mandelstam variables for massive kinematics with m3 = m4.
2114 : // tHG = tHat - m_gluino^2; uHG = uHat - m_gluino^2.
2115 : // (Note: tHG and uHG defined with opposite sign wrt Pythia 6)
2116 0 : tHG = -0.5 * (sH - tH + uH);
2117 0 : uHG = -0.5 * (sH + tH - uH);
2118 0 : tHG2 = tHG * tHG;
2119 0 : uHG2 = uHG * uHG;
2120 0 : s34Avg = 0.5 * (s3 + s4) - 0.25 * pow2(s3 - s4) / sH;
2121 :
2122 : // s-channel gluon contribution (only used if id1 == -id2)
2123 : // = Qss/s^2 in <Fuk11> including 2N*(N^2-1)/N^2 color factor.
2124 0 : sigS = 16./3. * (tHG2 + uHG2 + 2. * s34Avg * sH) / sH2;
2125 :
2126 0 : }
2127 :
2128 : //--------------------------------------------------------------------------
2129 :
2130 :
2131 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2132 :
2133 : double Sigma2qqbar2gluinogluino::sigmaHat() {
2134 :
2135 : // Only allow quark-antiquark incoming states
2136 0 : if (id1 * id2 > 0) return 0.0;
2137 :
2138 : // In-pair must both be up-type or both down-type
2139 0 : if ((id1+id2) % 2 != 0) return 0.0;
2140 :
2141 : // Flavor indices for the incoming quarks
2142 0 : int iQA = (abs(id1)+1)/2;
2143 0 : int iQB = (abs(id2)+1)/2;
2144 :
2145 : // Use up- or down-type squark-quark-gluino couplings
2146 0 : complex LsqqG[7][4];
2147 0 : complex RsqqG[7][4];
2148 0 : for (int iSq=1; iSq<=6; ++iSq) {
2149 0 : for (int iQ=1; iQ<=3; ++iQ) {
2150 0 : if (abs(id1) % 2 == 1) {
2151 0 : LsqqG[iSq][iQ] = coupSUSYPtr->LsddG[iSq][iQ];
2152 0 : RsqqG[iSq][iQ] = coupSUSYPtr->RsddG[iSq][iQ];
2153 0 : }
2154 : else {
2155 0 : LsqqG[iSq][iQ] = coupSUSYPtr->LsuuG[iSq][iQ];
2156 0 : RsqqG[iSq][iQ] = coupSUSYPtr->RsuuG[iSq][iQ];
2157 : }
2158 : }
2159 : }
2160 :
2161 : // sigHel contains (-1, 1), (1,-1), (-1,-1), ( 1, 1)
2162 : // divided by 4 for helicity average
2163 0 : vector<double> sigHel;
2164 0 : for (int iHel=0; iHel <4; ++iHel) sigHel.push_back(0.);
2165 :
2166 : // S-channel gluon contributions: only if id1 == -id2 (so iQA == iQB)
2167 0 : if (abs(id1) == abs(id2)) {
2168 : // S-channel squared
2169 0 : sigHel[0] += sigS;
2170 0 : sigHel[1] += sigS;
2171 0 : }
2172 :
2173 : // T, U, and S/T/U interferences
2174 0 : for (int iSq=1; iSq<=6; ++iSq) {
2175 0 : int idSq = ((iSq+2)/3)*1000000 + 2*((iSq-1)%3) + abs(id1-1) % 2 + 1;
2176 0 : double mSq2 = pow2(particleDataPtr->m0(idSq));
2177 : // Modified Mandelstam variables for massive kinematics with m3 = m4.
2178 : // tHG = tHat - m_gluino^2; uHG = uHat - m_gluino^2.
2179 0 : double tHsq = tHG + s34Avg - mSq2;
2180 0 : double uHsq = uHG + s34Avg - mSq2;
2181 :
2182 : // ST and SU interferences: only if id1 == -id2 (so iQA == iQB)
2183 : // incl 2N*(N^2 - 1)/N^2 color factor (note: original reference
2184 : // <Fuk11> was missing a factor 2 on the color factor here.)
2185 0 : if ( abs(id1) == abs(id2) ) {
2186 0 : double Qst1 = 16./3. * norm(LsqqG[iSq][iQA]) * (s34Avg * sH + tHG2);
2187 0 : double Qsu1 = 16./3. * norm(LsqqG[iSq][iQA]) * (s34Avg * sH + uHG2);
2188 0 : double Qst2 = 16./3. * norm(RsqqG[iSq][iQA]) * (s34Avg * sH + tHG2);
2189 0 : double Qsu2 = 16./3. * norm(RsqqG[iSq][iQA]) * (s34Avg * sH + uHG2);
2190 0 : double sigL = (Qst1 / tHsq + Qsu1 / uHsq) / sH;
2191 0 : double sigR = (Qst2 / tHsq + Qsu2 / uHsq) / sH;
2192 0 : sigHel[0] += sigL;
2193 0 : sigHel[1] += sigR;
2194 0 : }
2195 :
2196 : // T, U, and TU interferences
2197 0 : for (int jSq=1; jSq<=6; ++jSq) {
2198 0 : int idSqJ = ((jSq+2)/3)*1000000 + 2*((jSq-1)%3) + abs(id1-1) % 2 + 1;
2199 0 : double mSqJ2 = pow2(particleDataPtr->m0(idSqJ));
2200 : // Modified Mandelstam variables for massive kinematics with m3 = m4.
2201 : // tHG = tHat - m_gluino^2; uHG = uHat - m_gluino^2.
2202 0 : double tHsqJ = tHG + s34Avg - mSqJ2;
2203 0 : double uHsqJ = uHG + s34Avg - mSqJ2;
2204 :
2205 0 : double Q11 = real(LsqqG[iSq][iQA] * conj(LsqqG[iSq][iQB])
2206 0 : * conj(LsqqG[jSq][iQA]) * LsqqG[jSq][iQB]);
2207 0 : double Q12 = real(LsqqG[iSq][iQA] * conj(RsqqG[iSq][iQB])
2208 0 : * conj(LsqqG[jSq][iQA]) * RsqqG[jSq][iQB]);
2209 0 : double Q21 = real(RsqqG[iSq][iQA] * conj(LsqqG[iSq][iQB])
2210 0 : * conj(RsqqG[jSq][iQA]) * LsqqG[jSq][iQB]);
2211 0 : double Q22 = real(RsqqG[iSq][iQA] * conj(RsqqG[iSq][iQB])
2212 0 : * conj(RsqqG[jSq][iQA]) * RsqqG[jSq][iQB]);
2213 :
2214 0 : double Qtt11 = 64./27. * Q11 * tHG2;
2215 0 : double Qtt12 = 64./27. * Q12 * tHG2;
2216 0 : double Qtt21 = 64./27. * Q21 * tHG2;
2217 0 : double Qtt22 = 64./27. * Q22 * tHG2;
2218 :
2219 0 : double Quu11 = 64./27. * Q11 * uHG2;
2220 0 : double Quu12 = 64./27. * Q12 * uHG2;
2221 0 : double Quu21 = 64./27. * Q21 * uHG2;
2222 0 : double Quu22 = 64./27. * Q22 * uHG2;
2223 :
2224 0 : double Qtu11 = 16./27. * Q11 * (s34Avg * sH);
2225 0 : double Qtu12 = 16./27. * Q12 * (s34Avg * sH - tHG * uHG);
2226 0 : double Qtu21 = 16./27. * Q21 * (s34Avg * sH - tHG * uHG);
2227 0 : double Qtu22 = 16./27. * Q22 * (s34Avg * sH);
2228 :
2229 : // Cross sections for each helicity configuration (incl average fac 1/4)
2230 0 : sigHel[0] += Qtt11 / tHsq / tHsqJ
2231 0 : + Quu11 / uHsq / uHsqJ
2232 0 : + Qtu11 / tHsq / uHsqJ;
2233 0 : sigHel[1] += Qtt22 / tHsq / tHsqJ
2234 0 : + Quu22 / uHsq / uHsqJ
2235 0 : + Qtu22 / tHsq / uHsqJ;
2236 0 : sigHel[2] += Qtt12 / tHsq / tHsqJ
2237 0 : + Quu12 / uHsq / uHsqJ
2238 0 : + Qtu12 / tHsq / uHsqJ;
2239 0 : sigHel[3] += Qtt21 / tHsq / tHsqJ
2240 0 : + Quu21 / uHsq / uHsqJ
2241 0 : + Qtu21 / tHsq / uHsqJ;
2242 :
2243 : }
2244 :
2245 : }
2246 :
2247 : // Sum helicity contributions
2248 0 : double sigSum = sigHel[0] + sigHel[1] + sigHel[2] + sigHel[3];
2249 :
2250 : // Return 0 if all terms vanish, else compute and return cross section
2251 0 : if ( sigSum <= 0. ) return 0.0;
2252 :
2253 : // Answer
2254 0 : double sigma = (M_PI / 8. / sH2) * pow2(alpS) * sigSum * openFracPair;
2255 : return sigma;
2256 :
2257 0 : }
2258 :
2259 : //--------------------------------------------------------------------------
2260 :
2261 : // Select identity, colour and anticolour.
2262 :
2263 : void Sigma2qqbar2gluinogluino::setIdColAcol() {
2264 :
2265 : // Flavours are trivial.
2266 0 : setId( id1, id2, 1000021, 1000021);
2267 :
2268 : // Two colour flow topologies. Swap if first is antiquark.
2269 0 : if (rndmPtr->flat() < 0.5) setColAcol( 1, 0, 0, 2, 1, 3, 3, 2);
2270 0 : else setColAcol( 1, 0, 0, 2, 3, 2, 1, 3);
2271 0 : if (id1 < 0) swapColAcol();
2272 :
2273 0 : }
2274 :
2275 : //==========================================================================
2276 :
2277 : // Sigma1qq2antisquark
2278 : // R-parity violating squark production
2279 :
2280 : //--------------------------------------------------------------------------
2281 :
2282 : // Initialise process
2283 :
2284 : void Sigma1qq2antisquark::initProc(){
2285 :
2286 : //Typecast to the correct couplings
2287 0 : coupSUSYPtr = (CoupSUSY*) couplingsPtr;
2288 :
2289 : //Construct name of the process from lambda'' couplings
2290 :
2291 0 : nameSave = "q q' -> " + particleDataPtr->name(-idRes)+" + c.c";
2292 0 : codeSave = 2000 + 10*abs(idRes)/1000000 + abs(idRes)%10;
2293 0 : }
2294 :
2295 : //--------------------------------------------------------------------------
2296 :
2297 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2298 :
2299 : void Sigma1qq2antisquark::sigmaKin() {
2300 :
2301 : // Check if at least one RPV coupling non-zero
2302 0 : if(!coupSUSYPtr->isUDD) {
2303 0 : sigBW = 0.0;
2304 0 : return;
2305 : }
2306 :
2307 0 : mRes = particleDataPtr->m0(abs(idRes));
2308 0 : GammaRes = particleDataPtr->mWidth(abs(idRes));
2309 0 : m2Res = pow2(mRes);
2310 :
2311 0 : sigBW = sH * GammaRes/ ( pow2(sH - m2Res) + pow2(mRes * GammaRes) );
2312 0 : sigBW *= 2.0/3.0/mRes;
2313 :
2314 : // Width out only includes open channels.
2315 0 : widthOut = GammaRes * particleDataPtr->resOpenFrac(id3);
2316 0 : }
2317 :
2318 : //--------------------------------------------------------------------------
2319 :
2320 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2321 :
2322 : double Sigma1qq2antisquark::sigmaHat() {
2323 :
2324 : // Only allow (anti)quark-(anti)quark incoming states
2325 0 : if (id1*id2 <= 0) return 0.0;
2326 :
2327 : // Generation indices
2328 0 : int iA = (abs(id1)+1)/2;
2329 0 : int iB = (abs(id2)+1)/2;
2330 :
2331 : //Covert from pdg-code to ~u_i/~d_i basis
2332 0 : bool idown = (abs(idRes)%2 == 1) ? true : false;
2333 0 : int iC = (abs(idRes)/1000000 == 2)
2334 0 : ? (abs(idRes)%10+1)/2 + 3: (abs(idRes)%10+1)/2;
2335 :
2336 : // UDD structure
2337 0 : if (abs(id1)%2 == 0 && abs(id2)%2 == 0) return 0.0;
2338 0 : if (abs(id1)%2 == 1 && abs(id2)%2 == 1 && idown) return 0.0;
2339 0 : if ((abs(id1) + abs(id2))%2 == 1 && !idown) return 0.0;
2340 :
2341 : double sigma = 0.0;
2342 :
2343 0 : if(!idown){
2344 : // d_i d_j --> ~u*_k
2345 0 : for(int isq = 1; isq <=3; isq++){
2346 : // Loop over R-type squark contributions
2347 0 : sigma += pow2(coupSUSYPtr->rvUDD[isq][iA][iB])
2348 0 : * norm(coupSUSYPtr->Rusq[iC][isq+3]);
2349 : }
2350 0 : }else{
2351 : // u_i d_j --> d*_k
2352 : // Pick the correct coupling for d-u in-state
2353 0 : if(abs(id1)%2==1){
2354 0 : iA = (abs(id2)+1)/2;
2355 0 : iB = (abs(id1)+1)/2;
2356 0 : }
2357 0 : for(int isq = 1; isq <= 3; isq++){
2358 : // Loop over R-type squark contributions
2359 0 : sigma += pow2(coupSUSYPtr->rvUDD[iA][iB][isq])
2360 0 : * norm(coupSUSYPtr->Rdsq[iC][isq+3]);
2361 : }
2362 : }
2363 :
2364 0 : sigma *= sigBW;
2365 : return sigma;
2366 :
2367 0 : }
2368 :
2369 : //--------------------------------------------------------------------------
2370 :
2371 : // Select identity, colour and anticolour.
2372 :
2373 : void Sigma1qq2antisquark::setIdColAcol() {
2374 :
2375 : // Set flavours.
2376 0 : if(id1 < 0 && id2 < 0 ) setId( id1, id2, idRes);
2377 0 : else setId( id1, id2, -idRes);
2378 :
2379 : // Colour flow topologies. Swap when antiquarks.
2380 0 : if (abs(id1) < 9) setColAcol( 1, 0, 2, 0, 0, 3);
2381 0 : else setColAcol( 0, 0, 0, 0, 0, 0);
2382 0 : if (id1 < 0) swapColAcol();
2383 :
2384 0 : }
2385 :
2386 :
2387 : //==========================================================================
2388 :
2389 :
2390 : // Sigma2qqbar2chi0gluino
2391 : // Cross section for gaugino pair production: neutralino-gluino
2392 :
2393 : //--------------------------------------------------------------------------
2394 :
2395 : // Initialize process.
2396 :
2397 : void Sigma2qqbar2chi0gluino::initProc() {
2398 :
2399 : //Typecast to the correct couplings
2400 0 : coupSUSYPtr = (CoupSUSY*) couplingsPtr;
2401 :
2402 : // Construct name of process.
2403 0 : nameSave = "q qbar' -> " + particleDataPtr->name(id3) + " "
2404 0 : + particleDataPtr->name(id4);
2405 :
2406 : // Secondary open width fraction.
2407 0 : openFracPair = particleDataPtr->resOpenFrac(id3, id4);
2408 :
2409 0 : }
2410 :
2411 : //--------------------------------------------------------------------------
2412 :
2413 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2414 :
2415 : void Sigma2qqbar2chi0gluino::sigmaKin() {
2416 :
2417 : // Common flavour-independent factor.
2418 0 : sigma0 = M_PI * 4.0 / 9.0/ sH2 / coupSUSYPtr->sin2W * alpEM * alpS
2419 0 : * openFracPair;
2420 :
2421 : // Auxiliary factors for use below
2422 0 : ui = uH - s3;
2423 0 : uj = uH - s4;
2424 0 : ti = tH - s3;
2425 0 : tj = tH - s4;
2426 0 : }
2427 :
2428 : //--------------------------------------------------------------------------
2429 :
2430 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2431 :
2432 : double Sigma2qqbar2chi0gluino::sigmaHat() {
2433 :
2434 : // Only allow quark-antiquark incoming states
2435 0 : if (id1*id2 >= 0) return 0.0;
2436 :
2437 : // In-pair must both be up-type or both down-type
2438 0 : if ((id1+id2) % 2 != 0) return 0.0;
2439 :
2440 : // Swap T and U if antiquark-quark instead of quark-antiquark
2441 0 : if (id1<0) swapTU = true;
2442 :
2443 : // Shorthands
2444 0 : int idAbs1 = abs(id1);
2445 0 : int idAbs2 = abs(id2);
2446 :
2447 : // Flavour-dependent kinematics-dependent couplings.
2448 0 : complex QuLL(0.0),QtLL(0.0),QuRR(0.0),QtRR(0.0);
2449 0 : complex QuLR(0.0),QtLR(0.0),QuRL(0.0),QtRL(0.0);
2450 :
2451 : // Flavour indices
2452 0 : int ifl1 = (idAbs1+1) / 2;
2453 0 : int ifl2 = (idAbs2+1) / 2;
2454 :
2455 : complex (*LsddXloc)[4][6];
2456 : complex (*RsddXloc)[4][6];
2457 : complex (*LsuuXloc)[4][6];
2458 : complex (*RsuuXloc)[4][6];
2459 0 : LsddXloc = coupSUSYPtr->LsddX;
2460 0 : RsddXloc = coupSUSYPtr->RsddX;
2461 0 : LsuuXloc = coupSUSYPtr->LsuuX;
2462 0 : RsuuXloc = coupSUSYPtr->RsuuX;
2463 :
2464 : // Add t-channel squark flavour sums to QmXY couplings
2465 0 : for (int ksq=1; ksq<=6; ksq++) {
2466 :
2467 : // squark id and squark-subtracted u and t
2468 :
2469 : int idsq;
2470 0 : idsq=((ksq+2)/3)*1000000 + 2*((ksq-1) % 3) + (idAbs1+1) % 2 + 1;
2471 :
2472 0 : double msq2 = pow(particleDataPtr->m0(idsq),2);
2473 0 : double usq = uH - msq2;
2474 0 : double tsq = tH - msq2;
2475 :
2476 0 : complex Lsqq1X4;
2477 0 : complex Lsqq2X4;
2478 0 : complex Rsqq1X4;
2479 0 : complex Rsqq2X4;
2480 :
2481 0 : complex Lsqq1G;
2482 0 : complex Rsqq1G;
2483 0 : complex Lsqq2G;
2484 0 : complex Rsqq2G;
2485 :
2486 : // Couplings
2487 0 : Lsqq1X4 = LsuuXloc[ksq][ifl1][id4chi];
2488 0 : Lsqq2X4 = LsuuXloc[ksq][ifl2][id4chi];
2489 0 : Rsqq1X4 = RsuuXloc[ksq][ifl1][id4chi];
2490 0 : Rsqq2X4 = RsuuXloc[ksq][ifl2][id4chi];
2491 :
2492 0 : Lsqq1G = coupSUSYPtr->LsuuG[ksq][ifl1];
2493 0 : Rsqq1G = coupSUSYPtr->RsuuG[ksq][ifl1];
2494 0 : Lsqq2G = coupSUSYPtr->LsuuG[ksq][ifl2];
2495 0 : Rsqq2G = coupSUSYPtr->RsuuG[ksq][ifl2];
2496 :
2497 0 : if (idAbs1 % 2 != 0) {
2498 0 : Lsqq1X4 = LsddXloc[ksq][ifl1][id4chi];
2499 0 : Lsqq2X4 = LsddXloc[ksq][ifl2][id4chi];
2500 0 : Rsqq1X4 = RsddXloc[ksq][ifl1][id4chi];
2501 0 : Rsqq2X4 = RsddXloc[ksq][ifl2][id4chi];
2502 :
2503 0 : Lsqq1G = coupSUSYPtr->LsddG[ksq][ifl1];
2504 0 : Rsqq1G = coupSUSYPtr->RsddG[ksq][ifl1];
2505 0 : Lsqq2G = coupSUSYPtr->LsddG[ksq][ifl2];
2506 0 : Rsqq2G = coupSUSYPtr->RsddG[ksq][ifl2];
2507 0 : }
2508 :
2509 : // QuXY
2510 0 : QuLL += conj(Lsqq1X4)*Lsqq2G/usq;
2511 0 : QuRR += conj(Rsqq1X4)*Rsqq2G/usq;
2512 0 : QuLR += conj(Lsqq1X4)*Rsqq2G/usq;
2513 0 : QuRL += conj(Rsqq1X4)*Lsqq2G/usq;
2514 :
2515 : // QtXY
2516 0 : QtLL -= conj(Lsqq1G)*Lsqq2X4/tsq;
2517 0 : QtRR -= conj(Rsqq1G)*Rsqq2X4/tsq;
2518 0 : QtLR += conj(Lsqq1G)*Rsqq2X4/tsq;
2519 0 : QtRL += conj(Rsqq1G)*Lsqq2X4/tsq;
2520 :
2521 0 : }
2522 :
2523 : // Overall factor multiplying coupling
2524 0 : double fac = (1.0-coupSUSYPtr->sin2W);
2525 :
2526 : // Compute matrix element weight
2527 : double weight = 0;
2528 0 : double facLR = uH*tH - s3*s4;
2529 0 : double facMS = m3*m4*sH;
2530 :
2531 : // Average over separate helicity contributions
2532 : // LL (ha = -1, hb = +1) (divided by 4 for average)
2533 0 : weight += norm(QuLL) * ui * uj + norm(QtLL) * ti * tj
2534 0 : + 2 * real(conj(QuLL) * QtLL) * facMS;
2535 : // RR (ha = 1, hb = -1) (divided by 4 for average)
2536 0 : weight += norm(QtRR) * ti * tj + norm(QuRR) * ui * uj
2537 0 : + 2 * real(conj(QuRR) * QtRR) * facMS;
2538 : // RL (ha = 1, hb = 1) (divided by 4 for average)
2539 0 : weight += norm(QuRL) * ui * uj + norm(QtRL) * ti * tj
2540 0 : + real(conj(QuRL) * QtRL) * facLR;
2541 : // LR (ha = -1, hb = -1) (divided by 4 for average)
2542 0 : weight += norm(QuLR) * ui * uj + norm(QtLR) * ti * tj
2543 0 : + real(conj(QuLR) * QtLR) * facLR;
2544 :
2545 : // Cross section, including colour factor.
2546 0 : double sigma = sigma0 * weight / fac;
2547 :
2548 : // Answer.
2549 : return sigma;
2550 :
2551 0 : }
2552 :
2553 : //--------------------------------------------------------------------------
2554 :
2555 : // Select identity, colour and anticolour.
2556 :
2557 : void Sigma2qqbar2chi0gluino::setIdColAcol() {
2558 :
2559 : // Set flavours.
2560 0 : setId( id1, id2, id3, id4);
2561 :
2562 : // Colour flow topologies. Swap when antiquarks.
2563 0 : setColAcol( 1, 0, 0, 2, 1, 2, 0, 0);
2564 0 : if (id1 < 0) swapColAcol();
2565 :
2566 0 : }
2567 :
2568 : //==========================================================================
2569 :
2570 : // Sigma2qqbar2chargluino
2571 : // Cross section for gaugino pair production: chargino-gluino
2572 :
2573 : //--------------------------------------------------------------------------
2574 :
2575 : // Initialize process.
2576 :
2577 : void Sigma2qqbar2chargluino::initProc() {
2578 :
2579 : //Typecast to the correct couplings
2580 0 : coupSUSYPtr = (CoupSUSY*) couplingsPtr;
2581 :
2582 : // Construct name of process.
2583 0 : nameSave = "q qbar' -> " + particleDataPtr->name(id3) + " "
2584 0 : + particleDataPtr->name(id4) + " + c.c";
2585 :
2586 : // Secondary open width fraction.
2587 0 : openFracPair = particleDataPtr->resOpenFrac(id3, id4);
2588 :
2589 0 : }
2590 :
2591 : //--------------------------------------------------------------------------
2592 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2593 :
2594 : void Sigma2qqbar2chargluino::sigmaKin() {
2595 :
2596 : // Common flavour-independent factor.
2597 :
2598 0 : sigma0 = M_PI / sH2 * 4.0 / 9.0 / coupSUSYPtr->sin2W * alpEM * alpS ;
2599 0 : sigma0 /= 2.0 * (1 - coupSUSYPtr->sin2W) ;
2600 :
2601 : // Auxiliary factors for use below
2602 0 : ui = uH - s3;
2603 0 : uj = uH - s4;
2604 0 : ti = tH - s3;
2605 0 : tj = tH - s4;
2606 0 : }
2607 :
2608 : //--------------------------------------------------------------------------
2609 :
2610 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2611 :
2612 : double Sigma2qqbar2chargluino::sigmaHat() {
2613 :
2614 : // Only allow particle-antiparticle incoming states
2615 0 : if (id1*id2 >= 0) return 0.0;
2616 :
2617 : // Only allow incoming states with sum(charge) = final state
2618 0 : if (abs(id1) % 2 == abs(id2) % 2) return 0.0;
2619 0 : int isPos = (id4chi > 0 ? 1 : 0);
2620 0 : if (id1 < 0 && id1 > -19 && abs(id1) % 2 == 1-isPos ) return 0.0;
2621 0 : else if (id1 > 0 && id1 < 19 && abs(id1) % 2 == isPos ) return 0.0;
2622 :
2623 : // Flavour-dependent kinematics-dependent couplings.
2624 0 : int idAbs1 = abs(id1);
2625 0 : int iChar = abs(id4chi);
2626 :
2627 0 : complex QuLL(0.0),QtLL(0.0),QuRR(0.0),QtRR(0.0);
2628 0 : complex QuLR(0.0),QtLR(0.0),QuRL(0.0),QtRL(0.0);
2629 :
2630 : // Calculate everything from udbar -> ~chi+ ~chi0 template process
2631 0 : complex LsddGl;
2632 0 : complex RsddGl;
2633 0 : complex LsuuGl;
2634 0 : complex RsuuGl;
2635 : complex (*LsduXloc)[4][3];
2636 : complex (*RsduXloc)[4][3];
2637 : complex (*LsudXloc)[4][3];
2638 : complex (*RsudXloc)[4][3];
2639 :
2640 0 : LsduXloc = coupSUSYPtr->LsduX;
2641 0 : RsduXloc = coupSUSYPtr->RsduX;
2642 0 : LsudXloc = coupSUSYPtr->LsudX;
2643 0 : RsudXloc = coupSUSYPtr->RsudX;
2644 :
2645 : // u dbar , ubar d : do nothing
2646 : // dbar u , d ubar : swap 1<->2 and t<->u
2647 0 : int iGu = abs(id1)/2;
2648 0 : int iGd = (abs(id2)+1)/2;
2649 0 : if (idAbs1 % 2 != 0) {
2650 0 : swapTU = true;
2651 0 : iGu = abs(id2)/2;
2652 0 : iGd = (abs(id1)+1)/2;
2653 0 : }
2654 :
2655 : // Add t-channel squark flavour sums to QmXY couplings
2656 0 : for (int jsq=1; jsq<=6; jsq++) {
2657 :
2658 0 : int idsu=((jsq+2)/3)*1000000 + 2*((jsq-1) % 3) + 2 ;
2659 0 : int idsd=((jsq+2)/3)*1000000 + 2*((jsq-1) % 3) + 1 ;
2660 :
2661 0 : LsddGl = coupSUSYPtr->LsddG[jsq][iGd];
2662 0 : RsddGl = coupSUSYPtr->RsddG[jsq][iGd];
2663 0 : LsuuGl = coupSUSYPtr->LsuuG[jsq][iGu];
2664 0 : RsuuGl = coupSUSYPtr->RsuuG[jsq][iGu];
2665 :
2666 0 : double msd2 = pow(particleDataPtr->m0(idsd),2);
2667 0 : double msu2 = pow(particleDataPtr->m0(idsu),2);
2668 0 : double tsq = tH - msd2;
2669 0 : double usq = uH - msu2;
2670 :
2671 0 : QuLL += conj(LsuuGl) * conj(LsudXloc[jsq][iGd][iChar])/usq;
2672 0 : QuLR += conj(LsuuGl) * conj(RsudXloc[jsq][iGd][iChar])/usq;
2673 0 : QuRR += conj(RsuuGl) * conj(RsudXloc[jsq][iGd][iChar])/usq;
2674 0 : QuRL += conj(RsuuGl) * conj(LsudXloc[jsq][iGd][iChar])/usq;
2675 :
2676 0 : QtLL -= conj(LsduXloc[jsq][iGu][iChar]) * LsddGl/tsq;
2677 0 : QtRR -= conj(RsduXloc[jsq][iGu][iChar]) * RsddGl/tsq;
2678 0 : QtLR += conj(LsduXloc[jsq][iGu][iChar]) * RsddGl/tsq;
2679 0 : QtRL += conj(RsduXloc[jsq][iGu][iChar]) * LsddGl/tsq;
2680 0 : }
2681 :
2682 : // Compute matrix element weight
2683 : double weight = 0;
2684 :
2685 : // Average over separate helicity contributions
2686 : // (if swapped, swap ha, hb if computing polarized cross sections)
2687 : // LL (ha = -1, hb = +1) (divided by 4 for average)
2688 0 : weight += norm(QuLL) * ui * uj + norm(QtLL) * ti * tj
2689 0 : + 2 * real(conj(QuLL) * QtLL) * m3 * m4 * sH;
2690 : // RR (ha = 1, hb = -1) (divided by 4 for average)
2691 0 : weight += norm(QtRR) * ti * tj + norm(QuRR) * ui * uj
2692 0 : + 2 * real(conj(QuRR) * QtRR) * m3 * m4 * sH;
2693 : // RL (ha = 1, hb = 1) (divided by 4 for average)
2694 0 : weight += norm(QuRL) * ui * uj + norm(QtRL) * ti * tj
2695 0 : + real(conj(QuRL) * QtRL) * (uH * tH - s3 * s4);
2696 : // LR (ha = -1, hb = -1) (divided by 4 for average)
2697 0 : weight += norm(QuLR) * ui * uj + norm(QtLR) * ti * tj
2698 0 : + real(conj(QuLR) * QtLR) * (uH * tH - s3 * s4);
2699 :
2700 : // Cross section, including colour factor.
2701 0 : double sigma = sigma0 * weight;
2702 :
2703 : // Answer.
2704 : return sigma;
2705 :
2706 0 : }
2707 :
2708 : //--------------------------------------------------------------------------
2709 :
2710 : void Sigma2qqbar2chargluino::setIdColAcol() {
2711 :
2712 : // Set flavours.
2713 0 : setId( id1, id2, id3, id4);
2714 :
2715 : // Colour flow topologies. Swap when antiquarks.
2716 0 : setColAcol( 1, 0, 0, 2, 1, 2, 0, 0);
2717 0 : if (id1 < 0) swapColAcol();
2718 :
2719 0 : }
2720 :
2721 : //==========================================================================
2722 :
2723 : // Sigma2qqbar2sleptonantislepton
2724 : // Cross section for qqbar-initiated slepton-antislepton production
2725 :
2726 : //--------------------------------------------------------------------------
2727 :
2728 : // Initialize process.
2729 :
2730 : void Sigma2qqbar2sleptonantislepton::initProc() {
2731 :
2732 : //Typecast to the correct couplings
2733 0 : coupSUSYPtr = (CoupSUSY*) couplingsPtr;
2734 :
2735 : // Is this a ~e_i ~nu*_j, ~nu_i ~e*_j final state or ~e_i ~e*_j, ~nu_i ~nu*_j
2736 0 : if (abs(id3Sav) % 2 == abs(id4Sav) % 2) isUD = false;
2737 0 : else isUD = true;
2738 :
2739 : // Derive name
2740 0 : nameSave = "q qbar' -> "+particleDataPtr->name(abs(id3Sav))+" "+
2741 0 : particleDataPtr->name(-abs(id4Sav));
2742 0 : if (isUD) nameSave +=" + c.c.";
2743 :
2744 : // Extract isospin and mass-ordering indices
2745 :
2746 0 : if(isUD && abs(id3Sav)%2 == 0) {
2747 : // Make sure iGen3 is always slepton and iGen4 is always sneutrino
2748 0 : iGen3 = 3*(abs(id4Sav)/2000000) + (abs(id4Sav)%10+1)/2;
2749 0 : iGen4 = 3*(abs(id3Sav)/2000000) + (abs(id3Sav)%10+1)/2;
2750 0 : }
2751 : else {
2752 0 : iGen3 = 3*(abs(id3Sav)/2000000) + (abs(id3Sav)%10+1)/2;
2753 0 : iGen4 = 3*(abs(id4Sav)/2000000) + (abs(id4Sav)%10+1)/2;
2754 : }
2755 :
2756 : // Count 5 neutralinos in NMSSM
2757 0 : nNeut = (coupSUSYPtr->isNMSSM ? 5 : 4);
2758 :
2759 : // Store mass squares of all possible internal propagator lines;
2760 : // retained for future extension to leptonic initial states
2761 0 : m2Neut.resize(nNeut+1);
2762 0 : for (int iNeut=1;iNeut<=nNeut;iNeut++)
2763 0 : m2Neut[iNeut] = pow2(particleDataPtr->m0(coupSUSYPtr->idNeut(iNeut)));
2764 :
2765 : // Set sizes of some arrays to be used below
2766 0 : tNeut.resize(nNeut+1);
2767 0 : uNeut.resize(nNeut+1);
2768 :
2769 : // Shorthand for Weak mixing
2770 0 : xW = coupSUSYPtr->sin2W;
2771 :
2772 : // Secondary open width fraction.
2773 0 : openFracPair = particleDataPtr->resOpenFrac(id3Sav, id4Sav);
2774 :
2775 0 : }
2776 :
2777 : //--------------------------------------------------------------------------
2778 :
2779 : // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2780 :
2781 : void Sigma2qqbar2sleptonantislepton::sigmaKin() {
2782 :
2783 : // Z/W propagator
2784 0 : if (! isUD) {
2785 0 : double sV= sH - pow2(coupSUSYPtr->mZpole);
2786 0 : double d = pow2(sV) + pow2(coupSUSYPtr->mZpole * coupSUSYPtr->wZpole);
2787 0 : propZW = complex( sV / d, coupSUSYPtr->mZpole * coupSUSYPtr->wZpole / d);
2788 0 : } else {
2789 0 : double sV= sH - pow2(coupSUSYPtr->mWpole);
2790 0 : double d = pow2(sV) + pow2(coupSUSYPtr->mWpole * coupSUSYPtr->wWpole);
2791 0 : propZW = complex( sV / d, coupSUSYPtr->mWpole * coupSUSYPtr->wWpole / d);
2792 0 : }
2793 :
2794 : // Flavor-independent pre-factors
2795 0 : double comFacHat = M_PI/sH2 * openFracPair;
2796 :
2797 0 : sigmaEW = comFacHat * pow2(alpEM);
2798 0 : }
2799 :
2800 : //--------------------------------------------------------------------------
2801 :
2802 : // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2803 :
2804 : double Sigma2qqbar2sleptonantislepton::sigmaHat() {
2805 :
2806 : // In-pair must be opposite-sign
2807 0 : if (id1 * id2 > 0) return 0.0;
2808 :
2809 : // Check correct charge sum
2810 0 : if (isUD && abs(id1) %2 == abs(id2) % 2) return 0.0;
2811 0 : if (!isUD && abs(id1) % 2 != abs(id2) % 2) return 0.0;
2812 :
2813 : // No RH sneutrinos
2814 0 : if ( (abs(id3)%2 == 0 && abs(id3) > 2000000)
2815 0 : || (abs(id4)%2 == 0 && abs(id4) > 2000000) ) return 0.0;
2816 :
2817 : // Coded UD sigma is for udbar -> ~v~l'*. Swap t<->u for dbaru -> ~l~v*.
2818 0 : swapTU = (isUD && abs(id1) % 2 != 0);
2819 :
2820 : // Coded QQ sigma is for qqbar -> ~l~l*. Swap t<->u for qbarq -> ~l~l*.
2821 0 : if (!isUD && id1 < 0) swapTU = true;
2822 :
2823 : // Generation indices of incoming particles
2824 0 : int idIn1A = (swapTU) ? abs(id2) : abs(id1);
2825 0 : int idIn2A = (swapTU) ? abs(id1) : abs(id2);
2826 0 : int iGen1 = (idIn1A+1)/2;
2827 0 : int iGen2 = (idIn2A+1)/2;
2828 :
2829 : // Auxiliary factors for use below
2830 0 : for (int i=1; i<= nNeut; i++) {
2831 0 : tNeut[i] = tH - m2Neut[i];
2832 0 : uNeut[i] = uH - m2Neut[i];
2833 : }
2834 :
2835 0 : double eQ = (idIn1A % 2 == 0) ? 2./3. : -1./3. ;
2836 0 : double eSl = (abs(id3Sav) % 2 == 0) ? 0. : -1. ;
2837 :
2838 : // Initial values for pieces used for color-flow selection below
2839 0 : sumColS = 0.0;
2840 0 : sumColT = 0.0;
2841 0 : sumInterference = 0.0;
2842 :
2843 : // Common factor for LR and RL contributions
2844 0 : double facTU = uH*tH-s3*s4;
2845 :
2846 : // Opposite-isospin: udbar -> ~l~v*
2847 0 : if ( isUD ) {
2848 :
2849 : // s-channel W contribution (only contributes to LL helicities)
2850 0 : sumColS = sigmaEW / 32.0 / pow2(xW) / pow2(1.0-xW)
2851 0 : * norm(conj(coupSUSYPtr->LudW[iGen1][iGen2])
2852 0 : * coupSUSYPtr->LslsvW[iGen3][iGen4]) * facTU * norm(propZW);
2853 0 : }
2854 :
2855 : double CslZ;
2856 :
2857 : // s-channel Z/photon and interference
2858 0 : if (abs(id1) == abs(id2)) {
2859 :
2860 0 : CslZ = real(coupSUSYPtr->LslslZ[iGen3][iGen4]
2861 0 : + coupSUSYPtr->RslslZ[iGen3][iGen4]);
2862 0 : if (abs(id3)%2 == 0)
2863 0 : CslZ = real(coupSUSYPtr->LsvsvZ[iGen3][iGen4]
2864 0 : + coupSUSYPtr->RsvsvZ[iGen3][iGen4]);
2865 :
2866 : // gamma
2867 : // Factor 2 since contributes to both ha != hb helicities
2868 0 : sumColS += (abs(CslZ) > 0.0) ? 2. * pow2(eQ) * pow2(eSl) * sigmaEW
2869 0 : * facTU / pow2(sH) : 0.0;
2870 :
2871 : // Z/gamma interference
2872 0 : sumInterference += eQ * eSl * sigmaEW * facTU / 2.0 / xW / (1.-xW)
2873 0 : * sqrt(norm(propZW)) / sH * CslZ
2874 0 : * (coupSUSYPtr->LqqZ[idIn1A] + coupSUSYPtr->RqqZ[idIn1A]);
2875 :
2876 : // s-channel Z
2877 :
2878 0 : CslZ = norm(coupSUSYPtr->LslslZ[iGen3][iGen4]
2879 0 : + coupSUSYPtr->RslslZ[iGen3][iGen4]);
2880 0 : if (abs(id3Sav)%2 == 0)
2881 0 : CslZ = norm(coupSUSYPtr->LsvsvZ[iGen3][iGen4]
2882 0 : + coupSUSYPtr->RsvsvZ[iGen3][iGen4]);
2883 :
2884 0 : sumColS += sigmaEW * facTU / 16.0 / pow2(xW) / pow2(1.0-xW)
2885 0 : * norm(propZW) * CslZ
2886 0 : * ( pow2(coupSUSYPtr->LqqZ[idIn1A]) + pow2(coupSUSYPtr->RqqZ[idIn1A]) );
2887 0 : }
2888 :
2889 : // Cross section
2890 0 : double sigma = sumColS + sumColT + sumInterference;
2891 :
2892 : // Colour average
2893 0 : if(abs(id1) < 10) sigma /= 3.0;
2894 :
2895 : // Add cc term
2896 0 : if(isUD) sigma *= 2.0;
2897 :
2898 : // Return answer.
2899 : return sigma;
2900 :
2901 0 : }
2902 :
2903 : //--------------------------------------------------------------------------
2904 :
2905 : // Select identity, colour and anticolour.
2906 :
2907 : void Sigma2qqbar2sleptonantislepton::setIdColAcol() {
2908 :
2909 : // Set flavours.
2910 : int iSl, iSv;
2911 0 : if( isUD ){
2912 0 : iSl = (abs(id3)%2 == 0) ? abs(id3) : abs(id4);
2913 0 : iSv = (abs(id3)%2 == 0) ? abs(id4) : abs(id3);
2914 0 : if ((id1%2 + id2%2 ) > 0)
2915 0 : setId( id1, id2, -iSl, iSv);
2916 : else
2917 0 : setId( id1, id2, iSl, -iSv);
2918 : }
2919 : else
2920 0 : setId( id1, id2, abs(id3), -abs(id4));
2921 :
2922 0 : setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
2923 0 : if (id1 < 0 ) swapColAcol();
2924 :
2925 0 : }
2926 :
2927 : //==========================================================================
2928 :
2929 : } // end namespace Pythia8
|
