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