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