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) 2014 Torbjorn Sjostrand.
3 // Main authors of this file: N. Desai, P. Skands
4 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
5 // Please respect the MCnet Guidelines, see GUIDELINES for details.
6 
7 // Function definitions (not found in the header) for the
8 // supersymmetry simulation classes.
9 
10 #include "Pythia8/SigmaSUSY.h"
11 
12 namespace Pythia8 {
13 
14 //==========================================================================
15 
16 // Sigma2SUSY
17 // An intermediate class for SUSY 2 -> 2 with nontrivial decay angles.
18 
19 //--------------------------------------------------------------------------
20 
21 double Sigma2SUSY::weightDecay( Event& process, int iResBeg, int iResEnd) {
22 
23  // Do nothing if decays present already at input.
24 
25  // Identity of mother of decaying resonance(s).
26  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  sigma0 = M_PI / sH2 / coupSUSYPtr->sin2W * alpEM * alpS
712  * openFracPair;
713 
714  // Auxiliary factors for use below
715  ui = uH - s3;
716  uj = uH - s4;
717  ti = tH - s3;
718  tj = tH - s4;
719 
720 }
721 
722 //--------------------------------------------------------------------------
723 
724 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
725 
726 double Sigma2qg2chi0squark::sigmaHat() {
727 
728  // Antiquark -> antisquark
729  int idq = id1;
730  if (id1 == 21 || id1 == 22) idq = id2;
731  if (idq < 0) {
732  id4 = -abs(id4);
733  } else {
734  id4 = abs(id4);
735  }
736 
737  // tmp: only allow incoming quarks on side 1
738  // if (id1 < 0 || id1 == 21) return 0.0;
739 
740  // Generation index
741  int iGq = (abs(idq)+1)/2;
742 
743  // Only accept u(bar) -> ~u(bar) and d(bar) -> ~d(bar)
744  if (particleDataPtr->chargeType(idq) != particleDataPtr->chargeType(id4))
745  return 0.0;
746 
747  // Couplings
748  complex LsqqX, RsqqX;
749  if (idq % 2 == 0) {
750  LsqqX = coupSUSYPtr->LsuuX[id4sq][iGq][id3chi];
751  RsqqX = coupSUSYPtr->RsuuX[id4sq][iGq][id3chi];
752  }
753  else {
754  LsqqX = coupSUSYPtr->LsddX[id4sq][iGq][id3chi];
755  RsqqX = coupSUSYPtr->RsddX[id4sq][iGq][id3chi];
756  }
757 
758  // Prefactors : swap u and t if gq instead of qg
759  double fac1, fac2;
760  if (idq == id1) {
761  fac1 = -ui/sH + 2.0 * ( uH*tH - s4*s3 )/sH/tj;
762  fac2 = ti/tj * ( (tH + s4)/tj + (ti - uj)/sH );
763  } else {
764  fac1 = -ti/sH + 2.0 * ( uH*tH - s4*s3 )/sH/uj;
765  fac2 = ui/uj * ( (uH + s4)/uj + (ui - tj)/sH );
766  }
767 
768  // Compute matrix element weight
769  double weight = 0.0;
770 
771  // Average over separate helicity contributions
772  // (for qbar g : ha -> -ha )
773  // LL (ha = -1, hb = +1) (divided by 4 for average)
774  weight += fac2 * norm(LsqqX) / 2.0;
775  // RR (ha = 1, hb = -1) (divided by 4 for average)
776  weight += fac2 * norm(RsqqX) / 2.0;
777  // RL (ha = 1, hb = 1) (divided by 4 for average)
778  weight += fac2 * norm(RsqqX) / 2.0 + fac1 * norm(RsqqX);
779  // LR (ha = -1, hb = -1) (divided by 4 for average)
780  weight += fac2 * norm(LsqqX) / 2.0 + fac1 * norm(LsqqX);
781 
782  double sigma = sigma0 * weight;
783  if (abs(idq) < 9) sigma /= 3.;
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 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
838 
839 double Sigma2qg2charsquark::sigmaHat() {
840 
841  // Antiquark -> antisquark
842  int idq = (id1 == 21) ? id2 : id1;
843  if (idq > 0) {
844  id3 = id3Sav;
845  id4 = id4Sav;
846  } else {
847  id3 = -id3Sav;
848  id4 = -id4Sav;
849  }
850 
851  // Only accept u(bar) -> ~d(bar) and d(bar) -> ~u(bar)
852  if (particleDataPtr->chargeType(idq) == particleDataPtr->chargeType(id4))
853  return 0.0;
854 
855  // Generation index
856  int iGq = (abs(idq)+1)/2;
857 
858  // Couplings
859  complex LsqqX, RsqqX;
860  if (idq % 2 == 0) {
861  LsqqX = coupSUSYPtr->LsduX[id4sq][iGq][id3chi];
862  RsqqX = coupSUSYPtr->RsduX[id4sq][iGq][id3chi];
863  }
864  else {
865  LsqqX = coupSUSYPtr->LsduX[id4sq][iGq][id3chi];
866  RsqqX = coupSUSYPtr->RsduX[id4sq][iGq][id3chi];
867  }
868 
869  // Prefactors : swap u and t if gq instead of qg
870  double fac1, fac2;
871  if (idq == id1) {
872  fac1 = -ui/sH + 2.0 * ( uH*tH - s4*s3 )/sH/tj;
873  fac2 = ti/tj * ( (tH + s4)/tj + (ti - uj)/sH );
874  } else {
875  fac1 = -ti/sH + 2.0 * ( uH*tH - s4*s3 )/sH/uj;
876  fac2 = ui/uj * ( (uH + s4)/uj + (ui - tj)/sH );
877  }
878 
879  // Compute matrix element weight
880  double weight = 0.0;
881 
882  // Average over separate helicity contributions
883  // (a, b refers to qg configuration)
884  // LL (ha = -1, hb = +1) (divided by 4 for average)
885  weight += fac2 * norm(LsqqX) / 2.0;
886  // RR (ha = 1, hb = -1) (divided by 4 for average)
887  weight += fac2 * norm(RsqqX) / 2.0;
888  // RL (ha = 1, hb = 1) (divided by 4 for average)
889  weight += fac2 * norm(RsqqX) / 2.0 + fac1 * norm(RsqqX);
890  // LR (ha = -1, hb = -1) (divided by 4 for average)
891  weight += fac2 * norm(LsqqX) / 2.0 + fac1 * norm(LsqqX);
892 
893  double sigma = sigma0 * weight;
894  if (abs(idq) < 9) sigma /= 3.;
895 
896  // Answer.
897  return sigma * openFracPair;
898 
899 }
900 
901 //--------------------------------------------------------------------------
902 
903 // Select identity, colour and anticolour.
904 
905 void Sigma2qg2charsquark::setIdColAcol() {
906 
907  // Set flavours.
908  if (id1 > 0 && id2 > 0) {
909  setId( id1, id2, id3Sav, id4Sav);
910  } else {
911  setId( id1, id2,-id3Sav,-id4Sav);
912  }
913 
914  // Colour flow topology. Swap if first is gluon, or when antiquark.
915  if (id1 != 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
916  else setColAcol( 1, 2, 2, 0, 0, 0, 1, 0);
917  if (id1 < 0 || id2 < 0) swapColAcol();
918 
919 }
920 
921 //==========================================================================
922 
923 // Sigma2qq2squarksquark
924 // Cross section for squark-squark production
925 
926 //--------------------------------------------------------------------------
927 
928 // Initialize process.
929 
930 void Sigma2qq2squarksquark::initProc() {
931 
932  //Typecast to the correct couplings
933  coupSUSYPtr = (CoupSUSY*) couplingsPtr;
934 
935  // Extract mass-ordering indices
936  iGen3 = 3*(abs(id3Sav)/2000000) + (abs(id3Sav)%10+1)/2;
937  iGen4 = 3*(abs(id4Sav)/2000000) + (abs(id4Sav)%10+1)/2;
938 
939  // Is this a ~u_i ~d_j fial state or ~d_i ~d_j, ~u_i ~u_j
940  if (abs(id3Sav) % 2 == abs(id4Sav) % 2) isUD = false;
941  else isUD = true;
942 
943  // Derive name
944  nameSave = "q q' -> "+particleDataPtr->name(abs(id3Sav))+" "
945  +particleDataPtr->name(abs(id4Sav))+" + c.c.";
946 
947  // Count 5 neutralinos in NMSSM
948  nNeut = (coupSUSYPtr->isNMSSM ? 5 : 4);
949 
950  // Store mass squares of all possible internal propagator lines
951  m2Glu = pow2(particleDataPtr->m0(1000021));
952  m2Neut.resize(nNeut+1);
953  for (int iNeut=1;iNeut<=nNeut;iNeut++) {
954  m2Neut[iNeut] = pow2(particleDataPtr->m0(coupSUSYPtr->idNeut(iNeut)));
955  }
956  m2Char.resize(3);
957  m2Char[1] = pow2(particleDataPtr->m0(coupSUSYPtr->idChar(1)));
958  m2Char[2] = pow2(particleDataPtr->m0(coupSUSYPtr->idChar(2)));
959 
960  // Set sizes of some arrays to be used below
961  tNeut.resize(nNeut+1);
962  uNeut.resize(nNeut+1);
963  tChar.resize(3);
964  uChar.resize(3);
965 
966  // Secondary open width fraction.
967  openFracPair = particleDataPtr->resOpenFrac(id3Sav, id4Sav);
968 }
969 
970 //--------------------------------------------------------------------------
971 
972 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
973 
974 void Sigma2qq2squarksquark::sigmaKin() {
975 
976  // Weak mixing
977  double xW = coupSUSYPtr->sin2W;
978 
979  // pi/sH2
980  double comFacHat = M_PI/sH2 * openFracPair;
981 
982  // Channel-dependent but flavor-independent pre-factors
983  sigmaNeut = comFacHat * pow2(alpEM) / pow2(xW) / pow2(1-xW);
984  sigmaGlu = comFacHat * 2.0 * pow2(alpS) / 9.0;
985  if (isUD) {
986  sigmaChar = comFacHat * pow2(alpEM) / 4.0 / pow2(xW);
987  sigmaCharNeut = comFacHat * pow2(alpEM) / 3.0 / pow2(xW) / (1-xW);
988  sigmaCharGlu = comFacHat * 4.0 * alpEM * alpS / 9.0 / xW;
989  sigmaNeutGlu = 0.0;
990  } else {
991  sigmaChar = 0.0;
992  sigmaCharNeut = 0.0;
993  sigmaCharGlu = 0.0;
994  sigmaNeutGlu = comFacHat * 8.0 * alpEM * alpS / 9.0 / xW/(1-xW);
995  }
996 
997 }
998 
999 //--------------------------------------------------------------------------
1000 
1001 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1002 
1003 double Sigma2qq2squarksquark::sigmaHat() {
1004 
1005  // In-pair must be same-sign
1006  if (id1 * id2 < 0) return 0.0;
1007 
1008  // Check correct charge sum
1009  if (isUD && abs(id1) %2 == abs(id2) % 2) return 0.0;
1010  if (!isUD && abs(id1) % 2 != abs(id2) % 2) return 0.0;
1011  if (!isUD && abs(id1) % 2 != abs(id3Sav) % 2) return 0.0;
1012 
1013  // Coded sigma is for ud -> ~q~q'. Swap t and u for du -> ~q~q'.
1014  swapTU = (isUD && abs(id1) % 2 == 0);
1015  int idIn1A = (swapTU) ? abs(id2) : abs(id1);
1016  int idIn2A = (swapTU) ? abs(id1) : abs(id2);
1017 
1018  // Auxiliary factors for use below
1019  tGlu = tH - m2Glu;
1020  uGlu = uH - m2Glu;
1021  for (int i=1; i<= nNeut; i++) {
1022  tNeut[i] = tH - m2Neut[i];
1023  uNeut[i] = uH - m2Neut[i];
1024  if (isUD && i <= 2) {
1025  tChar[i] = tH - m2Char[i];
1026  uChar[i] = uH - m2Char[i];
1027  }
1028  }
1029 
1030  // Generation indices of incoming particles
1031  int iGen1 = (abs(idIn1A)+1)/2;
1032  int iGen2 = (abs(idIn2A)+1)/2;
1033 
1034  // Initial values for pieces used for color-flow selection below
1035  sumCt = 0.0;
1036  sumCu = 0.0;
1037  sumNt = 0.0;
1038  sumNu = 0.0;
1039  sumGt = 0.0;
1040  sumGu = 0.0;
1041  sumInterference = 0.0;
1042 
1043  // Common factor for LR and RL contributions
1044  double facTU = uH*tH-s3*s4;
1045 
1046  // Case A) Opposite-isospin: qq' -> ~d~u
1047  if ( isUD ) {
1048 
1049  // t-channel Charginos
1050  for (int k=1;k<=2;k++) {
1051 
1052  // Skip if only including gluinos
1053  if (settingsPtr->flag("SUSY:qq2squarksquark:onlyQCD")) continue;
1054 
1055  for (int l=1;l<=2;l++) {
1056 
1057  // kl-dependent factor for LL and RR contributions
1058  double facMS = sH*sqrt(m2Char[k]*m2Char[l]);
1059 
1060  // Note: Ckl defined as in [Boz07] with sigmaChar factored out
1061  // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
1062  complex Ckl[3][3];
1063  Ckl[1][1] = facMS
1064  * coupSUSYPtr->LsudX[iGen4][iGen2][k]
1065  * conj(coupSUSYPtr->LsduX[iGen3][iGen1][k])
1066  * conj(coupSUSYPtr->LsudX[iGen4][iGen2][l])
1067  * coupSUSYPtr->LsduX[iGen3][iGen1][l];
1068  Ckl[1][2] = facTU
1069  * coupSUSYPtr->RsudX[iGen4][iGen2][k]
1070  * conj(coupSUSYPtr->LsduX[iGen3][iGen1][k])
1071  * conj(coupSUSYPtr->RsudX[iGen4][iGen2][l])
1072  * coupSUSYPtr->LsduX[iGen3][iGen1][l];
1073  Ckl[2][1] = facTU
1074  * coupSUSYPtr->LsudX[iGen4][iGen2][k]
1075  * conj(coupSUSYPtr->RsduX[iGen3][iGen1][k])
1076  * conj(coupSUSYPtr->LsudX[iGen4][iGen2][l])
1077  * coupSUSYPtr->RsduX[iGen3][iGen1][l];
1078  Ckl[2][2] = facMS
1079  * coupSUSYPtr->RsudX[iGen4][iGen2][k]
1080  * conj(coupSUSYPtr->RsduX[iGen3][iGen1][k])
1081  * conj(coupSUSYPtr->RsudX[iGen4][iGen2][l])
1082  * coupSUSYPtr->RsduX[iGen3][iGen1][l];
1083 
1084  // Add to sum of t-channel charginos
1085  sumCt += sigmaChar * real(Ckl[1][1] + Ckl[1][2] + Ckl[2][1]
1086  + Ckl[2][2]) / tChar[k] / tChar[l];
1087 
1088  }
1089  }
1090 
1091  // u-channel Neutralinos
1092  for (int k=1;k<=nNeut;k++) {
1093 
1094  // Skip if only including gluinos
1095  if (settingsPtr->flag("SUSY:qq2squarksquark:onlyQCD")) continue;
1096 
1097  for (int l=1;l<=nNeut;l++) {
1098 
1099  // kl-dependent factor for LL, RR contributions
1100  double facMS = sH*sqrt(m2Neut[k]*m2Neut[l]);
1101 
1102  // Note: Nkl defined as in [Boz07] with sigmaNeut factored out
1103  // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
1104  complex Nkl[3][3];
1105  Nkl[1][1] = facMS
1106  * conj(coupSUSYPtr->LsuuX[iGen4][iGen1][k])
1107  * conj(coupSUSYPtr->LsddX[iGen3][iGen2][k])
1108  * coupSUSYPtr->LsuuX[iGen4][iGen1][l]
1109  * coupSUSYPtr->LsddX[iGen3][iGen2][l];
1110  Nkl[1][2] = facTU
1111  * conj(coupSUSYPtr->LsuuX[iGen4][iGen1][k])
1112  * conj(coupSUSYPtr->RsddX[iGen3][iGen2][k])
1113  * coupSUSYPtr->LsuuX[iGen4][iGen1][l]
1114  * coupSUSYPtr->RsddX[iGen3][iGen2][l];
1115  Nkl[2][1] = facTU
1116  * conj(coupSUSYPtr->RsuuX[iGen4][iGen1][k])
1117  * conj(coupSUSYPtr->LsddX[iGen3][iGen2][k])
1118  * coupSUSYPtr->RsuuX[iGen4][iGen1][l]
1119  * coupSUSYPtr->LsddX[iGen3][iGen2][l];
1120  Nkl[2][2] = facMS
1121  * conj(coupSUSYPtr->RsuuX[iGen4][iGen1][k])
1122  * conj(coupSUSYPtr->RsddX[iGen3][iGen2][k])
1123  * coupSUSYPtr->RsuuX[iGen4][iGen1][l]
1124  * coupSUSYPtr->RsddX[iGen3][iGen2][l];
1125 
1126  // Add to sum of u-channel neutralinos
1127  sumNu += sigmaNeut / uNeut[k] / uNeut[l]
1128  * real(Nkl[1][1] + Nkl[1][2] + Nkl[2][1] + Nkl[2][2]);
1129 
1130  }
1131  }
1132 
1133  // u-channel gluino
1134  // Note: Gij defined as in [Boz07] with sigmaGlu factored out
1135  // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
1136  double Gij[3][3];
1137  Gij[1][1] = norm(coupSUSYPtr->LsuuG[iGen4][iGen1]
1138  * coupSUSYPtr->LsddG[iGen3][iGen2]);
1139  Gij[1][2] = norm(coupSUSYPtr->LsuuG[iGen4][iGen1]
1140  * coupSUSYPtr->RsddG[iGen3][iGen2]);
1141  Gij[2][1] = norm(coupSUSYPtr->RsuuG[iGen4][iGen1]
1142  * coupSUSYPtr->LsddG[iGen3][iGen2]);
1143  Gij[2][2] = norm(coupSUSYPtr->RsuuG[iGen4][iGen1]
1144  * coupSUSYPtr->RsddG[iGen3][iGen2]);
1145  Gij[1][1] *= sH*m2Glu;
1146  Gij[1][2] *= facTU;
1147  Gij[2][1] *= facTU;
1148  Gij[2][2] *= sH*m2Glu;
1149  // Sum over polarizations
1150  sumGu += sigmaGlu * (Gij[1][1] + Gij[1][2] + Gij[2][1] + Gij[2][2])
1151  / pow2(uGlu);
1152 
1153  // chargino-neutralino interference
1154  for (int k=1;k<=2;k++) {
1155  for (int l=1;l<=nNeut;l++) {
1156 
1157  // Skip if only including gluinos
1158  if (settingsPtr->flag("SUSY:qq2squarksquark:onlyQCD")) continue;
1159 
1160  // Note: CNkl defined as in [Boz07] with pi/sH2 factored out
1161  // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
1162  double CNkl[3][3];
1163  CNkl[1][1] = real(coupSUSYPtr->LsudX[iGen4][iGen2][k]
1164  * conj(coupSUSYPtr->LsduX[iGen3][iGen1][k])
1165  * coupSUSYPtr->LsuuX[iGen4][iGen1][l]
1166  * coupSUSYPtr->LsddX[iGen3][iGen2][l]);
1167  CNkl[1][2] = real(coupSUSYPtr->RsudX[iGen4][iGen2][k]
1168  * conj(coupSUSYPtr->LsduX[iGen3][iGen1][k])
1169  * coupSUSYPtr->LsuuX[iGen4][iGen1][l]
1170  * coupSUSYPtr->RsddX[iGen3][iGen2][l]);
1171  CNkl[2][1] = real(coupSUSYPtr->LsudX[iGen4][iGen2][k]
1172  * conj(coupSUSYPtr->RsduX[iGen3][iGen1][k])
1173  * coupSUSYPtr->RsuuX[iGen4][iGen1][l]
1174  * coupSUSYPtr->LsddX[iGen3][iGen2][l]);
1175  CNkl[2][2] = real(coupSUSYPtr->RsudX[iGen4][iGen2][k]
1176  * conj(coupSUSYPtr->RsduX[iGen3][iGen1][k])
1177  * coupSUSYPtr->RsuuX[iGen4][iGen1][l]
1178  * coupSUSYPtr->RsddX[iGen3][iGen2][l]);
1179  CNkl[1][1] *= sH*sqrt(m2Char[k]*m2Neut[l]);
1180  CNkl[1][2] *= uH*tH-s3*s4;
1181  CNkl[2][1] *= uH*tH-s3*s4;
1182  CNkl[2][2] *= sH*sqrt(m2Char[k]*m2Neut[l]);
1183  // Sum over polarizations
1184  sumInterference += sigmaCharNeut * (CNkl[1][1] + CNkl[1][2]
1185  + CNkl[2][1] + CNkl[2][2]) / tChar[k] / uNeut[l];
1186  }
1187  }
1188 
1189  // chargino-gluino interference
1190  for (int k=1;k<=2;k++) {
1191 
1192  // Skip if only including gluinos
1193  if (settingsPtr->flag("SUSY:qq2squarksquark:onlyQCD")) continue;
1194 
1195  // Note: CGk defined as in [Boz07] with sigmaCharGlu factored out
1196  // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
1197  double CGk[3][3];
1198  CGk[1][1] = real(coupSUSYPtr->LsudX[iGen4][iGen2][k]
1199  * conj(coupSUSYPtr->LsduX[iGen3][iGen1][k])
1200  * conj(coupSUSYPtr->LsuuG[iGen4][iGen1])
1201  * conj(coupSUSYPtr->LsddG[iGen3][iGen2]));
1202  CGk[1][2] = real(coupSUSYPtr->RsudX[iGen4][iGen2][k]
1203  * conj(coupSUSYPtr->LsduX[iGen3][iGen1][k])
1204  * conj(coupSUSYPtr->LsuuG[iGen4][iGen1])
1205  * conj(coupSUSYPtr->RsddG[iGen3][iGen2]));
1206  CGk[2][1] = real(coupSUSYPtr->LsudX[iGen4][iGen2][k]
1207  * conj(coupSUSYPtr->RsduX[iGen3][iGen1][k])
1208  * conj(coupSUSYPtr->RsuuG[iGen4][iGen1])
1209  * conj(coupSUSYPtr->LsddG[iGen3][iGen2]));
1210  CGk[2][2] = real(coupSUSYPtr->RsudX[iGen4][iGen2][k]
1211  * conj(coupSUSYPtr->RsduX[iGen3][iGen1][k])
1212  * conj(coupSUSYPtr->RsuuG[iGen4][iGen1])
1213  * conj(coupSUSYPtr->RsddG[iGen3][iGen2]));
1214  CGk[1][1] *= sH*sqrt(m2Glu*m2Char[k]);
1215  CGk[1][2] *= uH*tH-s3*s4;
1216  CGk[2][1] *= uH*tH-s3*s4;
1217  CGk[2][2] *= sH*sqrt(m2Glu*m2Char[k]);
1218  // Sum over polarizations
1219  sumInterference += sigmaGlu * (CGk[1][1] + CGk[1][2] + CGk[2][1]
1220  + CGk[2][2]) / uGlu / tChar[k];
1221  }
1222  }
1223 
1224  // Case B) Same-isospin: qq' -> ~d~d , ~u~u
1225  else {
1226 
1227  // t-channel + u-channel Neutralinos + t/u interference
1228  for (int k=1;k<=nNeut;k++) {
1229 
1230  // Skip if only including gluinos
1231  if (settingsPtr->flag("SUSY:qq2squarksquark:onlyQCD")) continue;
1232 
1233  for (int l=1;l<=nNeut;l++) {
1234 
1235  // kl-dependent factor for LL and RR contributions
1236  double facMS = sH * particleDataPtr->m0(coupSUSYPtr->idNeut(k))
1237  * particleDataPtr->m0(coupSUSYPtr->idNeut(l));
1238 
1239  // Note: Nxkl defined as in [Boz07] with sigmaNeut factored out
1240  // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
1241  complex NTkl[3][3], NUkl[3][3], NTUkl[3][3];
1242  NTkl[1][1] = facMS
1243  * conj(coupSUSYPtr->getLsqqX(iGen4,idIn2A,k))
1244  * conj(coupSUSYPtr->getLsqqX(iGen3,idIn1A,k))
1245  * coupSUSYPtr->getLsqqX(iGen4,idIn2A,l)
1246  * coupSUSYPtr->getLsqqX(iGen3,idIn1A,l);
1247  NTkl[1][2] = facTU
1248  * conj(coupSUSYPtr->getRsqqX(iGen4,idIn2A,k))
1249  * conj(coupSUSYPtr->getLsqqX(iGen3,idIn1A,k))
1250  * coupSUSYPtr->getRsqqX(iGen4,idIn2A,l)
1251  * coupSUSYPtr->getLsqqX(iGen3,idIn1A,l);
1252  NTkl[2][1] = facTU
1253  * conj(coupSUSYPtr->getLsqqX(iGen4,idIn2A,k))
1254  * conj(coupSUSYPtr->getRsqqX(iGen3,idIn1A,k))
1255  * coupSUSYPtr->getLsqqX(iGen4,idIn2A,l)
1256  * coupSUSYPtr->getRsqqX(iGen3,idIn1A,l);
1257  NTkl[2][2] = facMS
1258  * conj(coupSUSYPtr->getRsqqX(iGen4,idIn2A,k))
1259  * conj(coupSUSYPtr->getRsqqX(iGen3,idIn1A,k))
1260  * coupSUSYPtr->getRsqqX(iGen4,idIn2A,l)
1261  * coupSUSYPtr->getRsqqX(iGen3,idIn1A,l);
1262  NUkl[1][1] = facMS
1263  * conj(coupSUSYPtr->getLsqqX(iGen3,idIn2A,k))
1264  * conj(coupSUSYPtr->getLsqqX(iGen4,idIn1A,k))
1265  * coupSUSYPtr->getLsqqX(iGen3,idIn2A,l)
1266  * coupSUSYPtr->getLsqqX(iGen4,idIn1A,l);
1267  NUkl[1][2] = facTU
1268  * conj(coupSUSYPtr->getRsqqX(iGen3,idIn2A,k))
1269  * conj(coupSUSYPtr->getLsqqX(iGen4,idIn1A,k))
1270  * coupSUSYPtr->getRsqqX(iGen3,idIn2A,l)
1271  * coupSUSYPtr->getLsqqX(iGen4,idIn1A,l);
1272  NUkl[2][1] = facTU
1273  * conj(coupSUSYPtr->getLsqqX(iGen3,idIn2A,k))
1274  * conj(coupSUSYPtr->getRsqqX(iGen4,idIn1A,k))
1275  * coupSUSYPtr->getLsqqX(iGen3,idIn2A,l)
1276  * coupSUSYPtr->getRsqqX(iGen4,idIn1A,l);
1277  NUkl[2][2] = facMS
1278  * conj(coupSUSYPtr->getRsqqX(iGen3,idIn2A,k))
1279  * conj(coupSUSYPtr->getRsqqX(iGen4,idIn1A,k))
1280  * coupSUSYPtr->getRsqqX(iGen3,idIn2A,l)
1281  * coupSUSYPtr->getRsqqX(iGen4,idIn1A,l);
1282  NTUkl[1][1] = facMS
1283  * real( conj(coupSUSYPtr->getLsqqX(iGen4,idIn2A,k))
1284  * conj(coupSUSYPtr->getLsqqX(iGen3,idIn1A,k))
1285  * coupSUSYPtr->getLsqqX(iGen3,idIn2A,l)
1286  * coupSUSYPtr->getLsqqX(iGen4,idIn1A,l) );
1287  NTUkl[1][2] = facTU
1288  * real( conj(coupSUSYPtr->getRsqqX(iGen4,idIn2A,k))
1289  * conj(coupSUSYPtr->getLsqqX(iGen3,idIn1A,k))
1290  * coupSUSYPtr->getRsqqX(iGen3,idIn2A,l)
1291  * coupSUSYPtr->getLsqqX(iGen4,idIn1A,l) );
1292  NTUkl[2][1] = facTU
1293  * real( conj(coupSUSYPtr->getLsqqX(iGen4,idIn2A,k))
1294  * conj(coupSUSYPtr->getRsqqX(iGen3,idIn1A,k))
1295  * coupSUSYPtr->getLsqqX(iGen3,idIn2A,l)
1296  * coupSUSYPtr->getRsqqX(iGen4,idIn1A,l) );
1297  NTUkl[2][2] = facMS
1298  * real( conj(coupSUSYPtr->getRsqqX(iGen4,idIn2A,k))
1299  * conj(coupSUSYPtr->getRsqqX(iGen3,idIn1A,k))
1300  * coupSUSYPtr->getRsqqX(iGen3,idIn2A,l)
1301  * coupSUSYPtr->getRsqqX(iGen4,idIn1A,l) );
1302 
1303  // Add to sums
1304  sumNt += sigmaNeut / tNeut[k] / tNeut[l]
1305  * real(NTkl[1][1] + NTkl[1][2] + NTkl[2][1] + NTkl[2][2]);
1306  sumNu += sigmaNeut / uNeut[k] / uNeut[l]
1307  * real(NUkl[1][1] + NUkl[1][2] + NUkl[2][1] + NUkl[2][2]);
1308  sumInterference += 2.0 / 3.0 * sigmaNeut
1309  * real(NTUkl[1][1] + NTUkl[1][2] + NTUkl[2][1] + NTUkl[2][2])
1310  / tNeut[k] / uNeut[l];
1311  }
1312 
1313  // Neutralino / Gluino interference
1314 
1315  // k-dependent factor for LL and RR contributions
1316  double facMS = sH * particleDataPtr->m0(coupSUSYPtr->idNeut(k))
1317  * particleDataPtr->m0(1000021);
1318 
1319  // Note: Nxkl defined as in [Boz07] with sigmaNeutGlu factored out
1320  // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
1321  complex NGA[3][3], NGB[3][3];
1322  NGA[1][1] = facMS
1323  * real( conj(coupSUSYPtr->getLsqqX(iGen4,idIn2A,k))
1324  * conj(coupSUSYPtr->getLsqqX(iGen3,idIn1A,k))
1325  * conj(coupSUSYPtr->getLsqqG(iGen3,idIn2A))
1326  * conj(coupSUSYPtr->getLsqqG(iGen4,idIn1A)) );
1327  NGA[1][2] = facTU
1328  * real( conj(coupSUSYPtr->getRsqqX(iGen4,idIn2A,k))
1329  * conj(coupSUSYPtr->getLsqqX(iGen3,idIn1A,k))
1330  * conj(coupSUSYPtr->getLsqqG(iGen3,idIn2A))
1331  * conj(coupSUSYPtr->getRsqqG(iGen4,idIn1A)) );
1332  NGA[2][1] = facTU
1333  * real( conj(coupSUSYPtr->getLsqqX(iGen4,idIn2A,k))
1334  * conj(coupSUSYPtr->getRsqqX(iGen3,idIn1A,k))
1335  * conj(coupSUSYPtr->getRsqqG(iGen3,idIn2A))
1336  * conj(coupSUSYPtr->getLsqqG(iGen4,idIn1A)) );
1337  NGA[2][2] = facMS
1338  * real( conj(coupSUSYPtr->getRsqqX(iGen4,idIn2A,k))
1339  * conj(coupSUSYPtr->getRsqqX(iGen3,idIn1A,k))
1340  * conj(coupSUSYPtr->getRsqqG(iGen3,idIn2A))
1341  * conj(coupSUSYPtr->getRsqqG(iGen4,idIn1A)) );
1342  NGB[1][1] = facMS
1343  * real( conj(coupSUSYPtr->getLsqqX(iGen3,idIn2A,k))
1344  * conj(coupSUSYPtr->getLsqqX(iGen4,idIn1A,k))
1345  * conj(coupSUSYPtr->getLsqqG(iGen4,idIn2A))
1346  * conj(coupSUSYPtr->getLsqqG(iGen3,idIn1A)) );
1347  NGB[1][2] = facMS
1348  * real( conj(coupSUSYPtr->getRsqqX(iGen3,idIn2A,k))
1349  * conj(coupSUSYPtr->getLsqqX(iGen4,idIn1A,k))
1350  * conj(coupSUSYPtr->getRsqqG(iGen4,idIn2A))
1351  * conj(coupSUSYPtr->getLsqqG(iGen3,idIn1A)) );
1352  NGB[2][1] = facMS
1353  * real( conj(coupSUSYPtr->getLsqqX(iGen3,idIn2A,k))
1354  * conj(coupSUSYPtr->getRsqqX(iGen4,idIn1A,k))
1355  * conj(coupSUSYPtr->getLsqqG(iGen4,idIn2A))
1356  * conj(coupSUSYPtr->getRsqqG(iGen3,idIn1A)) );
1357  NGB[2][2] = facMS
1358  * real( conj(coupSUSYPtr->getRsqqX(iGen3,idIn2A,k))
1359  * conj(coupSUSYPtr->getRsqqX(iGen4,idIn1A,k))
1360  * conj(coupSUSYPtr->getRsqqG(iGen4,idIn2A))
1361  * conj(coupSUSYPtr->getRsqqG(iGen3,idIn1A)) );
1362 
1363  // Add to sums
1364  sumInterference += sigmaNeutGlu *
1365  ( real(NGA[1][1] + NGA[1][2] + NGA[2][1] + NGA[2][2])
1366  / tNeut[k] / uGlu
1367  + real(NGB[1][1] + NGB[1][2] + NGB[2][1] + NGB[2][2])
1368  / uNeut[k] / tGlu );
1369  }
1370 
1371  // t-channel + u-channel Gluinos + t/u interference
1372 
1373  // factor for LL and RR contributions
1374  double facMS = sH * m2Glu;
1375 
1376  // Note: GT, GU defined as in [Boz07] with sigmaGlu factored out
1377  // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
1378  complex GT[3][3], GU[3][3], GTU[3][3];
1379  GT[1][1] = facMS
1380  * norm(coupSUSYPtr->getLsqqG(iGen4,idIn2A)
1381  * coupSUSYPtr->getLsqqG(iGen3,idIn1A));
1382  GT[1][2] = facTU
1383  * norm(coupSUSYPtr->getRsqqG(iGen4,idIn2A)
1384  * coupSUSYPtr->getLsqqG(iGen3,idIn1A));
1385  GT[2][1] = facTU
1386  * norm(coupSUSYPtr->getLsqqG(iGen4,idIn2A)
1387  * coupSUSYPtr->getRsqqG(iGen3,idIn1A));
1388  GT[2][2] = facMS
1389  * norm(coupSUSYPtr->getRsqqG(iGen4,idIn2A)
1390  * coupSUSYPtr->getRsqqG(iGen3,idIn1A));
1391  GU[1][1] = facMS
1392  * norm(coupSUSYPtr->getLsqqG(iGen3,idIn2A)
1393  * coupSUSYPtr->getLsqqG(iGen4,idIn1A));
1394  GU[1][2] = facTU
1395  * norm(coupSUSYPtr->getLsqqG(iGen3,idIn2A)
1396  * coupSUSYPtr->getRsqqG(iGen4,idIn1A));
1397  GU[2][1] = facTU
1398  * norm(coupSUSYPtr->getRsqqG(iGen3,idIn2A)
1399  * coupSUSYPtr->getLsqqG(iGen4,idIn1A));
1400  GU[2][2] = facMS
1401  * norm(coupSUSYPtr->getRsqqG(iGen3,idIn2A)
1402  * coupSUSYPtr->getRsqqG(iGen4,idIn1A));
1403 
1404  GTU[1][1] = facMS
1405  * real(coupSUSYPtr->getLsqqG(iGen3,idIn1A)
1406  * coupSUSYPtr->getLsqqG(iGen4,idIn2A)
1407  * conj(coupSUSYPtr->getLsqqG(iGen3,idIn2A))
1408  * conj(coupSUSYPtr->getLsqqG(iGen4,idIn1A)) );
1409 
1410  GTU[1][2] = facTU
1411  * real(coupSUSYPtr->getLsqqG(iGen3,idIn1A)
1412  * coupSUSYPtr->getRsqqG(iGen4,idIn2A)
1413  * conj(coupSUSYPtr->getRsqqG(iGen3,idIn2A))
1414  * conj(coupSUSYPtr->getLsqqG(iGen4,idIn1A)) );
1415 
1416  GTU[2][1] = facTU
1417  * real(coupSUSYPtr->getRsqqG(iGen3,idIn1A)
1418  * coupSUSYPtr->getLsqqG(iGen4,idIn2A)
1419  * conj(coupSUSYPtr->getLsqqG(iGen3,idIn2A))
1420  * conj(coupSUSYPtr->getRsqqG(iGen4,idIn1A)) );
1421 
1422  GTU[2][2] = facMS
1423  * real(coupSUSYPtr->getRsqqG(iGen3,idIn1A)
1424  * coupSUSYPtr->getRsqqG(iGen4,idIn2A)
1425  * conj(coupSUSYPtr->getRsqqG(iGen3,idIn2A))
1426  * conj(coupSUSYPtr->getRsqqG(iGen4,idIn1A)) );
1427 
1428  // Add to sums
1429  sumGt += sigmaGlu * real(GT[1][1] + GT[1][2] + GT[2][1] + GT[2][2])
1430  / pow2(tGlu) ;
1431  sumGu += sigmaGlu * real(GU[1][1] + GU[1][2] + GU[2][1] + GU[2][2])
1432  / pow2(uGlu) ;
1433  sumInterference += - 2.0 / 3.0 * sigmaGlu
1434  * real(GTU[1][1] + GTU[1][2] + GTU[2][1] + GTU[2][2])
1435  / tGlu / uGlu;
1436 
1437  }
1438 
1439  // Cross section
1440  double sigma = sumNt + sumNu + sumCt + sumCu + sumGt + sumGu
1441  + sumInterference;
1442 
1443  // Identical particles?
1444  if (id3Sav == id4Sav) sigma /= 2.0;
1445 
1446  // Return answer.
1447  return sigma;
1448 
1449 }
1450 
1451 //--------------------------------------------------------------------------
1452 
1453 // Select identity, colour and anticolour.
1454 
1455 void Sigma2qq2squarksquark::setIdColAcol() {
1456 
1457  // Set flavours.
1458  if (id1 > 0 && id2 > 0) {
1459  setId( id1, id2, id3Sav, id4Sav);
1460  } else {
1461  // 1,2 -> -3,-4
1462  setId( id1, id2,-id3Sav,-id4Sav);
1463  }
1464 
1465  // Coded sigma is for ud -> ~q~q'. Swap t and u for du -> ~q~q'.
1466  swapTU = (isUD && abs(id1) % 2 == 0);
1467 
1468  // Select colour flow topology
1469  // Recompute contributions to this particular in- out- flavour combination
1470  sigmaHat();
1471  // A: t-channel neutralino, t-channel chargino, or u-channel gluino
1472  double sumA = sumNt + sumCt + sumGu;
1473  double sumAB = sumNt + sumNu + sumCt + sumCu + sumGt + sumGu;
1474  if (swapTU) sumA = sumAB - sumA;
1475  setColAcol( 1, 0, 2, 0, 1, 0, 2, 0);
1476  // B: t-channel gluino or u-channel neutralino
1477  if (rndmPtr->flat()*sumAB > sumA) setColAcol( 1, 0, 2, 0, 2, 0, 1, 0);
1478 
1479  // Switch to anti-colors if antiquarks
1480  if (id1 < 0 || id2 < 0) swapColAcol();
1481 
1482 }
1483 
1484 //==========================================================================
1485 
1486 // Sigma2qqbar2squarkantisquark
1487 // Cross section for qqbar-initiated squark-antisquark production
1488 
1489 //--------------------------------------------------------------------------
1490 
1491 // Initialize process.
1492 
1493 void Sigma2qqbar2squarkantisquark::initProc() {
1494 
1495  //Typecast to the correct couplings
1496  coupSUSYPtr = (CoupSUSY*) couplingsPtr;
1497 
1498  // Is this a ~u_i ~d*_j, ~d_i ~u*_j final state or ~d_i ~d*_j, ~u_i ~u*_j
1499  if (abs(id3Sav) % 2 == abs(id4Sav) % 2) isUD = false;
1500  else isUD = true;
1501 
1502  // Extract isospin and mass-ordering indices
1503  if(isUD && abs(id3Sav)%2 == 1){
1504  iGen3 = 3*(abs(id4Sav)/2000000) + (abs(id3Sav)%10+1)/2;
1505  iGen4 = 3*(abs(id3Sav)/2000000) + (abs(id4Sav)%10+1)/2;
1506  }
1507  else {
1508  iGen3 = 3*(abs(id3Sav)/2000000) + (abs(id3Sav)%10+1)/2;
1509  iGen4 = 3*(abs(id4Sav)/2000000) + (abs(id4Sav)%10+1)/2;
1510  }
1511 
1512  // Derive name
1513  nameSave = "q qbar' -> "+particleDataPtr->name(abs(id3Sav))+" "+
1514  particleDataPtr->name(-abs(id4Sav));
1515  if (isUD && abs(id3Sav) != abs(id4Sav)) nameSave +=" + c.c.";
1516 
1517  // Count 5 neutralinos in NMSSM
1518  nNeut = (coupSUSYPtr->isNMSSM ? 5 : 4);
1519 
1520  // Store mass squares of all possible internal propagator lines
1521  m2Glu = pow2(particleDataPtr->m0(1000021));
1522  m2Neut.resize(nNeut+1);
1523  for (int iNeut=1;iNeut<=nNeut;iNeut++)
1524  m2Neut[iNeut] = pow2(particleDataPtr->m0(coupSUSYPtr->idNeut(iNeut)));
1525 
1526  // Set sizes of some arrays to be used below
1527  tNeut.resize(nNeut+1);
1528  uNeut.resize(nNeut+1);
1529 
1530  // Shorthand for Weak mixing
1531  xW = coupSUSYPtr->sin2W;
1532 
1533  // Secondary open width fraction.
1534  openFracPair = particleDataPtr->resOpenFrac(id3Sav, id4Sav);
1535 
1536 }
1537 
1538 //--------------------------------------------------------------------------
1539 
1540 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
1541 
1542 void Sigma2qqbar2squarkantisquark::sigmaKin() {
1543 
1544  // Z/W propagator
1545  if (! isUD) {
1546  double sV= sH - pow2(coupSUSYPtr->mZpole);
1547  double d = pow2(sV) + pow2(coupSUSYPtr->mZpole * coupSUSYPtr->wZpole);
1548  propZW = complex( sV / d, coupSUSYPtr->mZpole * coupSUSYPtr->wZpole / d);
1549  } else {
1550  double sV= sH - pow2(coupSUSYPtr->mWpole);
1551  double d = pow2(sV) + pow2(coupSUSYPtr->mWpole * coupSUSYPtr->wWpole);
1552  propZW = complex( sV / d, coupSUSYPtr->mWpole * coupSUSYPtr->wWpole / d);
1553  }
1554 
1555  // Flavor-independent pre-factors
1556  double comFacHat = M_PI/sH2 * openFracPair;
1557 
1558  sigmaEW = comFacHat * pow2(alpEM);
1559  sigmaGlu = comFacHat * 2.0 * pow2(alpS) / 9.0;
1560  sigmaEWG = comFacHat * 8.0 * alpEM * alpS / 9.0;
1561 
1562 }
1563 
1564 //--------------------------------------------------------------------------
1565 
1566 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
1567 
1568 double Sigma2qqbar2squarkantisquark::sigmaHat() {
1569 
1570  // In-pair must be opposite-sign
1571  if (id1 * id2 > 0) return 0.0;
1572 
1573  // Check correct charge sum
1574  if (isUD && abs(id1) %2 == abs(id2) % 2) return 0.0;
1575  if (!isUD && abs(id1) % 2 != abs(id2) % 2) return 0.0;
1576 
1577  // Check if using QCD diagrams only
1578  bool onlyQCD = settingsPtr->flag("SUSY:qqbar2squarkantisquark:onlyQCD");
1579 
1580  // Coded UD sigma is for udbar -> ~u~d'*. Swap t<->u for dbaru -> ~u~d'*.
1581  swapTU = (isUD && abs(id1) % 2 != 0);
1582 
1583  // Coded QQ sigma is for qqbar -> ~q~q*. Swap t<->u for qbarq -> ~q~q*.
1584  if (!isUD && id1 < 0) swapTU = true;
1585 
1586  // Generation indices of incoming particles
1587  int idIn1A = (swapTU) ? abs(id2) : abs(id1);
1588  int idIn2A = (swapTU) ? abs(id1) : abs(id2);
1589  int iGen1 = (idIn1A+1)/2;
1590  int iGen2 = (idIn2A+1)/2;
1591 
1592  // Auxiliary factors for use below
1593  tGlu = tH - m2Glu;
1594  uGlu = uH - m2Glu;
1595  for (int i=1; i<= nNeut; i++) {
1596  tNeut[i] = tH - m2Neut[i];
1597  uNeut[i] = uH - m2Neut[i];
1598  }
1599 
1600  // Initial values for pieces used for color-flow selection below
1601  sumColS = 0.0;
1602  sumColT = 0.0;
1603  sumInterference = 0.0;
1604 
1605  // Common factor for LR and RL contributions
1606  double facTU = uH*tH-s3*s4;
1607 
1608  // Case A) Opposite-isospin: udbar -> ~u~d*
1609  if ( isUD ) {
1610 
1611  // s-channel W contribution (only contributes to LL helicities)
1612  if ( !onlyQCD ) {
1613  sumColS += sigmaEW / 16.0 / pow2(xW) / pow2(1.0-xW)
1614  * norm(conj(coupSUSYPtr->LudW[iGen1][iGen2])
1615  * coupSUSYPtr->LsusdW[iGen3][iGen4]) * facTU
1616  * norm(propZW);
1617  }
1618 
1619  // t-channel gluino contributions
1620  double GT[3][3];
1621  double facLR = m2Glu * sH;
1622  // LL, LR, RL, RR
1623  GT[1][1] = facTU * norm(conj(coupSUSYPtr->LsddG[iGen4][iGen2])
1624  *coupSUSYPtr->LsuuG[iGen3][iGen1]);
1625  GT[1][2] = facLR * norm(conj(coupSUSYPtr->RsddG[iGen4][iGen2])
1626  *coupSUSYPtr->LsuuG[iGen3][iGen1]);
1627  GT[2][1] = facLR * norm(conj(coupSUSYPtr->LsddG[iGen4][iGen2])
1628  *coupSUSYPtr->RsuuG[iGen3][iGen1]);
1629  GT[2][2] = facTU * norm(conj(coupSUSYPtr->RsddG[iGen4][iGen2])
1630  *coupSUSYPtr->RsuuG[iGen3][iGen1]);
1631  // leading color flow for t-channel gluino is annihilation-like
1632  sumColS += sigmaGlu / pow2(tGlu)
1633  * (GT[1][1] + GT[1][2] + GT[2][1] + GT[2][2]);
1634 
1635  // W-Gluino interference (only contributes to LL helicities)
1636  if ( !onlyQCD ) {
1637  sumColS += sigmaEWG / 4.0 / xW / (1-xW)
1638  * real(conj(coupSUSYPtr->LsuuG[iGen3][iGen1])
1639  * coupSUSYPtr->LsddG[iGen4][iGen2]
1640  * conj(coupSUSYPtr->LudW[iGen1][iGen2])
1641  * coupSUSYPtr->LsusdW[iGen3][iGen4]) * facTU
1642  / tGlu * sqrt(norm(propZW));
1643  }
1644 
1645  // t-channel neutralinos
1646  // NOT YET IMPLEMENTED !
1647 
1648  }
1649 
1650  // Case B) Same-isospin: qqbar -> ~d~d* , ~u~u*
1651  else {
1652 
1653  double eQ = (idIn1A % 2 == 0) ? 2./3. : 1./3. ;
1654  double eSq = (abs(id3Sav) % 2 == 0) ? 2./3. : 1./3. ;
1655 
1656  // s-channel gluon (strictly flavor-diagonal)
1657  if (abs(id3Sav) == abs(id4Sav) && abs(id1) == abs(id2)) {
1658  // Factor 2 since contributes to both ha != hb helicities
1659  sumColT += 2. * sigmaGlu * facTU / pow2(sH);
1660  }
1661 
1662  // t-channel gluino (only for in-isospin = out-isospin).
1663  if (eQ == eSq) {
1664  // Sum over helicities.
1665  double GT[3][3];
1666  double facLR = sH * m2Glu;
1667  GT[1][1] = facTU * norm(coupSUSYPtr->getLsqqG(iGen3,idIn1A)
1668  * conj(coupSUSYPtr->getLsqqG(iGen4,idIn2A)));
1669  GT[1][2] = facLR * norm(coupSUSYPtr->getLsqqG(iGen3,idIn1A)
1670  * conj(coupSUSYPtr->getRsqqG(iGen4,idIn2A)));
1671  GT[2][1] = facLR * norm(coupSUSYPtr->getRsqqG(iGen3,idIn1A)
1672  * conj(coupSUSYPtr->getLsqqG(iGen4,idIn2A)));
1673  GT[2][2] = facTU * norm(coupSUSYPtr->getRsqqG(iGen3,idIn1A)
1674  * conj(coupSUSYPtr->getRsqqG(iGen4,idIn2A)));
1675  // Add contribution to color topology: S
1676  sumColS += sigmaGlu / pow2(tGlu)
1677  * ( GT[1][1] + GT[1][2] + GT[2][1] + GT[2][2]);
1678 
1679  // gluon-gluino interference (strictly flavor-diagonal)
1680  if (abs(id3Sav) == abs(id4Sav) && abs(id1) == abs(id2)) {
1681  double GG11, GG22;
1682  GG11 = - facTU * 2./3.
1683  * real( conj(coupSUSYPtr->getLsqqG(iGen3,idIn1A))
1684  * coupSUSYPtr->getLsqqG(iGen4,idIn2A));
1685  GG22 = - facTU * 2./3.
1686  * real( conj(coupSUSYPtr->getRsqqG(iGen3,idIn1A))
1687  * coupSUSYPtr->getRsqqG(iGen4,idIn2A));
1688  // Sum over two contributing helicities
1689  sumInterference += sigmaGlu / sH / tGlu
1690  * ( GG11 + GG22 );
1691  }
1692 
1693  }
1694 
1695  // Skip the rest if only including QCD diagrams
1696  if (onlyQCD) return sumColT+sumColS+sumInterference;
1697 
1698  // s-channel photon (strictly flavor-diagonal) and Z/gamma interference
1699  if (abs(id3Sav) == abs(id4Sav) && abs(id1) == abs(id2)) {
1700 
1701  // gamma
1702  // Factor 2 since contributes to both ha != hb helicities
1703  sumColS += 2. * pow2(eQ) * pow2(eSq) * sigmaEW * facTU / pow2(sH);
1704 
1705  // Z/gamma interference
1706  double CsqZ = real(coupSUSYPtr->LsusuZ[iGen3][iGen4]
1707  + coupSUSYPtr->RsusuZ[iGen3][iGen4]);
1708  if (abs(id3Sav)%2 != 0) CsqZ = real(coupSUSYPtr->LsdsdZ[iGen3][iGen4]
1709  + coupSUSYPtr->RsdsdZ[iGen3][iGen4]);
1710  sumColS += eQ * eSq * sigmaEW * facTU / 2.0 / xW / (1.-xW)
1711  * sqrt(norm(propZW)) / sH * CsqZ
1712  * (coupSUSYPtr->LqqZ[idIn1A] + coupSUSYPtr->LqqZ[idIn2A]);
1713 
1714  // Gluino/gamma interference (only for same-isospin)
1715  if (eQ == eSq) {
1716  double CsqG11 = real(conj(coupSUSYPtr->LsuuG[iGen3][iGen1])
1717  *coupSUSYPtr->LsuuG[iGen4][iGen2]);
1718  double CsqG22 = real(conj(coupSUSYPtr->RsuuG[iGen3][iGen1])
1719  *coupSUSYPtr->RsuuG[iGen4][iGen2]);
1720  if (id3Sav%2 != 0) {
1721  CsqG11 = real(conj(coupSUSYPtr->LsddG[iGen3][iGen1])
1722  *coupSUSYPtr->LsddG[iGen4][iGen2]);
1723  CsqG22 = real(conj(coupSUSYPtr->RsddG[iGen3][iGen1])
1724  *coupSUSYPtr->RsddG[iGen4][iGen2]);
1725  }
1726  sumColS += eQ * eSq * sigmaEWG * facTU
1727  * (CsqG11 + CsqG22) / sH / tGlu;
1728  }
1729  }
1730 
1731  // s-channel Z (only for q flavor = qbar flavor)
1732  if (abs(id1) == abs(id2)) {
1733  double CsqZ = norm(coupSUSYPtr->LsusuZ[iGen3][iGen4]
1734  + coupSUSYPtr->RsusuZ[iGen3][iGen4]);
1735  if (abs(id3Sav)%2 != 0) CsqZ = norm(coupSUSYPtr->LsdsdZ[iGen3][iGen4]
1736  + coupSUSYPtr->RsdsdZ[iGen3][iGen4]);
1737  sumColS += sigmaEW * facTU / 16.0 / pow2(xW) / pow2(1.0-xW)
1738  * norm(propZW) * CsqZ * ( pow2(coupSUSYPtr->LqqZ[idIn1A])
1739  + pow2(coupSUSYPtr->RqqZ[idIn1A]) );
1740 
1741  // Z/gluino interference (only for in-isospin = out-isospin)
1742  if (eQ == eSq) {
1743  double GZ11 = real(conj(coupSUSYPtr->getLsqqG(iGen3,idIn1A))
1744  *coupSUSYPtr->getLsqqG(iGen4,idIn2A)
1745  *(coupSUSYPtr->getLsqsqZ(id3Sav,id4Sav)
1746  +coupSUSYPtr->getRsqsqZ(id3Sav,id4Sav)))
1747  *coupSUSYPtr->LqqZ[idIn1A];
1748  double GZ22 = real(conj(coupSUSYPtr->getRsqqG(iGen3,idIn1A))
1749  *coupSUSYPtr->getRsqqG(iGen4,idIn2A)
1750  *(coupSUSYPtr->getLsqsqZ(id3Sav,id4Sav)
1751  +coupSUSYPtr->getRsqsqZ(id3Sav,id4Sav)))
1752  *coupSUSYPtr->RqqZ[idIn1A];
1753  sumColS += sigmaEWG * facTU / 4.0 / xW / (1.-xW)
1754  * ( GZ11 + GZ22 ) * sqrt(norm(propZW)) / tGlu;
1755  }
1756  }
1757 
1758  // t-channel neutralinos
1759  // NOT YET IMPLEMENTED !
1760 
1761  }
1762 
1763  // Cross section
1764  double sigma = sumColS + sumColT + sumInterference;
1765 
1766  // Return answer.
1767  return sigma;
1768 
1769 }
1770 
1771 //--------------------------------------------------------------------------
1772 
1773 // Select identity, colour and anticolour.
1774 
1775 void Sigma2qqbar2squarkantisquark::setIdColAcol() {
1776 
1777  // Check if charge conjugate final state?
1778  isCC = false;
1779  if (isUD && ( (id1-1)%2 < 0 || (id2-1)%2 < 0 )) isCC = true;
1780 
1781  //check if charge conjugate
1782  id3 = (isCC) ? -id3Sav : id3Sav;
1783  id4 = (isCC) ? -id4Sav : id4Sav;
1784 
1785  // Set flavours.
1786  setId( id1, id2, id3, id4);
1787 
1788  // Coded UD sigma is for udbar -> ~u~d'*. Swap t<->u for dbaru -> ~u~d'*.
1789  // Coded QQ sigma is for qqbar -> ~q~q*. Swap t<->u for qbarq -> ~q~q*.
1790  if (isUD) {
1791  swapTU = (abs(id1) % 2 != 0);
1792  } else {
1793  swapTU = (id1 < 0);
1794  }
1795 
1796  // Select colour flow topology
1797  // Recompute individual contributions to this in-out flavour combination
1798  sigmaHat();
1799  double R = rndmPtr->flat();
1800  double fracS = sumColS / (sumColS + sumColT) ;
1801  // S: color flow as in S-channel singlet
1802  if (R < fracS) {
1803  setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
1804  if (swapTU) setColAcol( 0, 1, 1, 0, 2, 0, 0, 2);
1805  }
1806  // T: color flow as in T-channel singlet
1807  else {
1808  setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
1809  if (swapTU) setColAcol( 0, 1, 2, 0, 2, 0, 0, 1);
1810  }
1811 
1812  if (isCC) swapColAcol();
1813 
1814 }
1815 
1816 //==========================================================================
1817 
1818 // Sigma2gg2squarkantisquark
1819 // Cross section for gg-initiated squark-antisquark production
1820 
1821 //--------------------------------------------------------------------------
1822 
1823 // Initialize process.
1824 
1825 void Sigma2gg2squarkantisquark::initProc() {
1826 
1827  //Typecast to the correct couplings
1828  coupSUSYPtr = (CoupSUSY*) couplingsPtr;
1829 
1830  // Process Name
1831  nameSave = "g g -> "+particleDataPtr->name(abs(id3Sav))+" "
1832  +particleDataPtr->name(-abs(id4Sav));
1833 
1834  // Squark pole mass
1835  m2Sq = pow2(particleDataPtr->m0(id3Sav));
1836 
1837  // Secondary open width fraction.
1838  openFracPair = particleDataPtr->resOpenFrac(id3Sav, id4Sav);
1839 
1840 }
1841 
1842 //--------------------------------------------------------------------------
1843 
1844 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
1845 
1846 void Sigma2gg2squarkantisquark::sigmaKin() {
1847 
1848  // Modified Mandelstam variables for massive kinematics with m3 = m4.
1849  // tHSq = tHat - m_squark^2; uHSq = uHat - m_squark^2.
1850  // double s34Avg = 0.5 * (s3 + s4) - 0.25 * pow2(s3 - s4) / sH;
1851  double tHSq = -0.5 * (sH - tH + uH);
1852  double uHSq = -0.5 * (sH + tH - uH);
1853  // ! (NEED TO CHECK THAT THESE APPLIED CORRECTLY BELOW) !
1854  // ! (PRELIMINARY CROSS-CHECKS WITH PYTHIA 6 COME OUT OK) !
1855 
1856  // Helicity-independent prefactor
1857  double comFacHat = M_PI/sH2 * pow2(alpS) / 128.0
1858  * ( 24.0 * (1.0 - 2*tHSq*uHSq/sH2) - 8.0/3.0 );
1859 
1860  // Helicity-dependent factors
1861  sigma = 0.0;
1862  for (int ha=-1;ha<=1;ha += 2) {
1863  for (int hb=-1;hb<=1;hb += 2) {
1864  // Divide by 4 for helicity average
1865  sigma += comFacHat / 4.0
1866  * ( (1.0-ha*hb)
1867  - 2.0 * sH*m2Sq/tHSq/uHSq
1868  * ( 1.0 - ha*hb - sH*m2Sq/tHSq/uHSq));
1869  }
1870  }
1871 
1872 }
1873 
1874 //--------------------------------------------------------------------------
1875 
1876 // Select identity, colour and anticolour.
1877 
1878 void Sigma2gg2squarkantisquark::setIdColAcol() {
1879 
1880  // Set flavours.
1881  setId( id1, id2, id3Sav, id4Sav);
1882 
1883  // Set color flow (random for now)
1884  double R = rndmPtr->flat();
1885  if (R < 0.5) setColAcol( 1, 2, 2, 3, 1, 0, 0, 3);
1886  else setColAcol( 1, 2, 3, 1, 3, 0, 0, 2);
1887 
1888 }
1889 
1890 //==========================================================================
1891 
1892 // Sigma2qg2squarkgluino
1893 // Cross section for squark-gluino production
1894 
1895 //--------------------------------------------------------------------------
1896 
1897 // Initialize process.
1898 
1899 void Sigma2qg2squarkgluino::initProc() {
1900 
1901  //Typecast to the correct couplings
1902  coupSUSYPtr = (CoupSUSY*) couplingsPtr;
1903 
1904  // Derive name
1905  nameSave = "q g -> "+particleDataPtr->name(abs(id3Sav))+" gluino + c.c.";
1906 
1907  // Final-state mass squares
1908  m2Glu = pow2(particleDataPtr->m0(1000021));
1909  m2Sq = pow2(particleDataPtr->m0(id3Sav));
1910 
1911  // Secondary open width fraction.
1912  openFracPair = particleDataPtr->resOpenFrac(id3Sav, 1000021);
1913 
1914 }
1915 
1916 //--------------------------------------------------------------------------
1917 
1918 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
1919 
1920 void Sigma2qg2squarkgluino::sigmaKin() {
1921 
1922  // Common pre-factor
1923  comFacHat = (M_PI / sH2) * pow2(alpS) * 0.5 * openFracPair;
1924 
1925  // Invariants (still with Pythia 6 sign convention)
1926  double tGlu = m2Glu-tH;
1927  double uGlu = m2Glu-uH;
1928  double tSq = m2Sq-tH;
1929  double uSq = m2Sq-uH;
1930 
1931  // Color flow A: quark color annihilates with anticolor of g
1932  sigmaA = 0.5*4./9.* tGlu/sH + (tGlu*sH+2.*m2Glu*tSq)/pow2(tGlu) -
1933  ( (sH-m2Sq+m2Glu)*(-tSq)-sH*m2Glu )/sH/(-tGlu)
1934  + 0.5*1./2.*( tSq*(tH+2.*uH+m2Glu)-tGlu*(sH-2.*tSq)
1935  + (-uGlu)*(tH+m2Glu+2.*m2Sq) )/2./tGlu/uSq;
1936  // Color flow B: quark and gluon colors iterchanged
1937  sigmaB = 4./9.*(-uGlu)*(uH+m2Sq)/pow2(uSq)
1938  + 1./18.* (sH*(uH+m2Glu) + 2.*(m2Sq-m2Glu)*uGlu)/sH/(-uSq)
1939  + 0.5*4./9.*tGlu/sH
1940  + 0.5*1./2.*(tSq*(tH+2.*uH+m2Glu)-tGlu*(sH-2.*tSq)
1941  + (-uGlu)*(tH+m2Glu+2.*m2Sq))/2./tGlu/uSq;
1942 
1943 }
1944 
1945 //--------------------------------------------------------------------------
1946 
1947 double Sigma2qg2squarkgluino::sigmaHat() {
1948 
1949  // Check whether right incoming flavor
1950  int idQA = (id1 == 21) ? id2 : id1;
1951  int idSq = (abs(id3) == 10000021) ? id4 : id3;
1952 
1953  // Check for charge conservation
1954  if(idQA%2 != idSq%2) return 0.0;
1955 
1956  int idQ = (abs(idQA)+1)/2;
1957  idSq = 3 * (id3Sav / 2000000) + (id3Sav % 10 + 1)/2;
1958 
1959  double mixingFac;
1960  if(abs(idQA) % 2 == 1)
1961  mixingFac = norm(coupSUSYPtr->LsddG[idSq][idQ])
1962  + norm(coupSUSYPtr->RsddG[idSq][idQ]);
1963  else
1964  mixingFac = norm(coupSUSYPtr->LsuuG[idSq][idQ])
1965  + norm(coupSUSYPtr->RsuuG[idSq][idQ]);
1966 
1967  return mixingFac * comFacHat * (sigmaA + sigmaB);
1968 }
1969 
1970 //--------------------------------------------------------------------------
1971 
1972 // Select identity, colour and anticolour.
1973 
1974 void Sigma2qg2squarkgluino::setIdColAcol() {
1975 
1976  // Check if charge conjugate final state?
1977  int idQ = (id1 == 21) ? id2 : id1;
1978  id3 = (idQ > 0) ? id3Sav : -id3Sav;
1979  id4 = 1000021;
1980 
1981  // Set flavors
1982  setId( id1, id2, id3, id4);
1983 
1984  // Select color flow A or B (see above)
1985  // Recompute individual contributions to this in-out flavour combination
1986  sigmaHat();
1987  double R = rndmPtr->flat()*(sigmaA+sigmaB);
1988  if (idQ == id1) {
1989  setColAcol(1,0,2,1,3,0,2,3);
1990  if (R > sigmaA) setColAcol(1,0,2,3,2,0,1,3);
1991  } else {
1992  setColAcol(2,1,1,0,3,0,2,3);
1993  if (R > sigmaB) setColAcol(2,3,1,0,2,0,1,3);
1994  }
1995  if (idQ < 0) swapColAcol();
1996 
1997  // Use reflected kinematics if gq initial state
1998  if (id1 == 21) swapTU = true;
1999 
2000 }
2001 
2002 //==========================================================================
2003 
2004 // Sigma2gg2gluinogluino
2005 // Cross section for gluino pair production from gg initial states
2006 // (validated against Pythia 6)
2007 
2008 //--------------------------------------------------------------------------
2009 
2010 // Initialize process.
2011 
2012 void Sigma2gg2gluinogluino::initProc() {
2013 
2014  //Typecast to the correct couplings
2015  coupSUSYPtr = (CoupSUSY*) couplingsPtr;
2016 
2017  // Secondary open width fraction.
2018  openFracPair = particleDataPtr->resOpenFrac(1000021, 1000021);
2019 
2020 }
2021 
2022 //--------------------------------------------------------------------------
2023 
2024 // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence.
2025 
2026 void Sigma2gg2gluinogluino::sigmaKin() {
2027 
2028  // Modified Mandelstam variables for massive kinematics with m3 = m4.
2029  // tHG = tHat - m_gluino^2; uHG = uHat - m_gluino^2.
2030  double s34Avg = 0.5 * (s3 + s4) - 0.25 * pow2(s3 - s4) / sH;
2031  double tHG = -0.5 * (sH - tH + uH);
2032  double uHG = -0.5 * (sH + tH - uH);
2033  double tHG2 = tHG * tHG;
2034  double uHG2 = uHG * uHG;
2035 
2036  // Calculate kinematics dependence.
2037  sigTS = (tHG * uHG - 2. * s34Avg * (tHG + 2. * s34Avg)) / tHG2
2038  + (tHG * uHG + s34Avg * (uHG - tHG)) / (sH * tHG);
2039  sigUS = (tHG * uHG - 2. * s34Avg * (uHG + 2. * s34Avg)) / uHG2
2040  + (tHG * uHG + s34Avg * (tHG - uHG)) / (sH * uHG);
2041  sigTU = 2. * tHG * uHG / sH2 + s34Avg * (sH - 4. * s34Avg)
2042  / (tHG * uHG);
2043  sigSum = sigTS + sigUS + sigTU;
2044 
2045  // Answer contains factor 1/2 from identical gluinos.
2046  sigma = (M_PI / sH2) * pow2(alpS) * (9./4.) * 0.5 * sigSum
2047  * openFracPair;
2048 
2049 }
2050 
2051 //--------------------------------------------------------------------------
2052 
2053 // Select identity, colour and anticolour.
2054 
2055 void Sigma2gg2gluinogluino::setIdColAcol() {
2056 
2057  // Flavours are trivial.
2058  setId( id1, id2, 1000021, 1000021);
2059 
2060  // Three colour flow topologies, each with two orientations.
2061  double sigRand = sigSum * rndmPtr->flat();
2062  if (sigRand < sigTS) setColAcol( 1, 2, 2, 3, 1, 4, 4, 3);
2063  else if (sigRand < sigTS + sigUS)
2064  setColAcol( 1, 2, 3, 1, 3, 4, 4, 2);
2065  else setColAcol( 1, 2, 3, 4, 1, 4, 3, 2);
2066  if (rndmPtr->flat() > 0.5) swapColAcol();
2067 
2068 }
2069 
2070 //==========================================================================
2071 
2072 // Sigma2qqbar2gluinogluino
2073 // Cross section for gluino pair production from qqbar initial states
2074 // (validated against Pythia 6 for SLHA1 case)
2075 
2076 //--------------------------------------------------------------------------
2077 
2078 // Initialize process.
2079 
2080 void Sigma2qqbar2gluinogluino::initProc() {
2081 
2082  //Typecast to the correct couplings
2083  coupSUSYPtr = (CoupSUSY*) couplingsPtr;
2084 
2085  // Secondary open width fraction.
2086  openFracPair = particleDataPtr->resOpenFrac(1000021, 1000021);
2087 
2088 }
2089 
2090 //--------------------------------------------------------------------------
2091 
2092 // Begin evaluate d(sigmaHat)/d(tHat); flavour-independent part.
2093 
2094 void Sigma2qqbar2gluinogluino::sigmaKin() {
2095 
2096  // Modified Mandelstam variables for massive kinematics with m3 = m4.
2097  // tHG = tHat - m_gluino^2; uHG = uHat - m_gluino^2.
2098  // (Note: tHG and uHG defined with opposite sign wrt Pythia 6)
2099  tHG = -0.5 * (sH - tH + uH);
2100  uHG = -0.5 * (sH + tH - uH);
2101  tHG2 = tHG * tHG;
2102  uHG2 = uHG * uHG;
2103  s34Avg = 0.5 * (s3 + s4) - 0.25 * pow2(s3 - s4) / sH;
2104 
2105  // s-channel gluon contribution (only used if id1 == -id2)
2106  // = Qss/s^2 in <Fuk11> including 2N*(N^2-1)/N^2 color factor.
2107  sigS = 16./3. * (tHG2 + uHG2 + 2. * s34Avg * sH) / sH2;
2108 
2109 }
2110 
2111 //--------------------------------------------------------------------------
2112 
2113 
2114 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2115 
2116 double Sigma2qqbar2gluinogluino::sigmaHat() {
2117 
2118  // Only allow quark-antiquark incoming states
2119  if (id1 * id2 > 0) return 0.0;
2120 
2121  // In-pair must both be up-type or both down-type
2122  if ((id1+id2) % 2 != 0) return 0.0;
2123 
2124  // Flavor indices for the incoming quarks
2125  int iQA = (abs(id1)+1)/2;
2126  int iQB = (abs(id2)+1)/2;
2127 
2128  // Use up- or down-type squark-quark-gluino couplings
2129  complex LsqqG[7][4];
2130  complex RsqqG[7][4];
2131  for (int iSq=1; iSq<=6; ++iSq) {
2132  for (int iQ=1; iQ<=3; ++iQ) {
2133  if (abs(id1) % 2 == 1) {
2134  LsqqG[iSq][iQ] = coupSUSYPtr->LsddG[iSq][iQ];
2135  RsqqG[iSq][iQ] = coupSUSYPtr->RsddG[iSq][iQ];
2136  }
2137  else {
2138  LsqqG[iSq][iQ] = coupSUSYPtr->LsuuG[iSq][iQ];
2139  RsqqG[iSq][iQ] = coupSUSYPtr->RsuuG[iSq][iQ];
2140  }
2141  }
2142  }
2143 
2144  // sigHel contains (-1, 1), (1,-1), (-1,-1), ( 1, 1)
2145  // divided by 4 for helicity average
2146  vector<double> sigHel;
2147  for (int iHel=0; iHel <4; ++iHel) sigHel.push_back(0.);
2148 
2149  // S-channel gluon contributions: only if id1 == -id2 (so iQA == iQB)
2150  if (abs(id1) == abs(id2)) {
2151  // S-channel squared
2152  sigHel[0] += sigS;
2153  sigHel[1] += sigS;
2154  }
2155 
2156  // T, U, and S/T/U interferences
2157  for (int iSq=1; iSq<=6; ++iSq) {
2158  int idSq = ((iSq+2)/3)*1000000 + 2*((iSq-1)%3) + abs(id1-1) % 2 + 1;
2159  double mSq2 = pow2(particleDataPtr->m0(idSq));
2160  // Modified Mandelstam variables for massive kinematics with m3 = m4.
2161  // tHG = tHat - m_gluino^2; uHG = uHat - m_gluino^2.
2162  double tHsq = tHG + s34Avg - mSq2;
2163  double uHsq = uHG + s34Avg - mSq2;
2164 
2165  // ST and SU interferences: only if id1 == -id2 (so iQA == iQB)
2166  // incl 2N*(N^2 - 1)/N^2 color factor (note: original reference
2167  // <Fuk11> was missing a factor 2 on the color factor here.)
2168  if ( abs(id1) == abs(id2) ) {
2169  double Qst1 = 16./3. * norm(LsqqG[iSq][iQA]) * (s34Avg * sH + tHG2);
2170  double Qsu1 = 16./3. * norm(LsqqG[iSq][iQA]) * (s34Avg * sH + uHG2);
2171  double Qst2 = 16./3. * norm(RsqqG[iSq][iQA]) * (s34Avg * sH + tHG2);
2172  double Qsu2 = 16./3. * norm(RsqqG[iSq][iQA]) * (s34Avg * sH + uHG2);
2173  double sigL = (Qst1 / tHsq + Qsu1 / uHsq) / sH;
2174  double sigR = (Qst2 / tHsq + Qsu2 / uHsq) / sH;
2175  sigHel[0] += sigL;
2176  sigHel[1] += sigR;
2177  }
2178 
2179  // T, U, and TU interferences
2180  for (int jSq=1; jSq<=6; ++jSq) {
2181  int idSqJ = ((jSq+2)/3)*1000000 + 2*((jSq-1)%3) + abs(id1-1) % 2 + 1;
2182  double mSqJ2 = pow2(particleDataPtr->m0(idSqJ));
2183  // Modified Mandelstam variables for massive kinematics with m3 = m4.
2184  // tHG = tHat - m_gluino^2; uHG = uHat - m_gluino^2.
2185  double tHsqJ = tHG + s34Avg - mSqJ2;
2186  double uHsqJ = uHG + s34Avg - mSqJ2;
2187 
2188  double Q11 = real(LsqqG[iSq][iQA] * conj(LsqqG[iSq][iQB])
2189  * conj(LsqqG[jSq][iQA]) * LsqqG[jSq][iQB]);
2190  double Q12 = real(LsqqG[iSq][iQA] * conj(RsqqG[iSq][iQB])
2191  * conj(LsqqG[jSq][iQA]) * RsqqG[jSq][iQB]);
2192  double Q21 = real(RsqqG[iSq][iQA] * conj(LsqqG[iSq][iQB])
2193  * conj(RsqqG[jSq][iQA]) * LsqqG[jSq][iQB]);
2194  double Q22 = real(RsqqG[iSq][iQA] * conj(RsqqG[iSq][iQB])
2195  * conj(RsqqG[jSq][iQA]) * RsqqG[jSq][iQB]);
2196 
2197  double Qtt11 = 64./27. * Q11 * tHG2;
2198  double Qtt12 = 64./27. * Q12 * tHG2;
2199  double Qtt21 = 64./27. * Q21 * tHG2;
2200  double Qtt22 = 64./27. * Q22 * tHG2;
2201 
2202  double Quu11 = 64./27. * Q11 * uHG2;
2203  double Quu12 = 64./27. * Q12 * uHG2;
2204  double Quu21 = 64./27. * Q21 * uHG2;
2205  double Quu22 = 64./27. * Q22 * uHG2;
2206 
2207  double Qtu11 = 16./27. * Q11 * (s34Avg * sH);
2208  double Qtu12 = 16./27. * Q12 * (s34Avg * sH - tHG * uHG);
2209  double Qtu21 = 16./27. * Q21 * (s34Avg * sH - tHG * uHG);
2210  double Qtu22 = 16./27. * Q22 * (s34Avg * sH);
2211 
2212  // Cross sections for each helicity configuration (incl average fac 1/4)
2213  sigHel[0] += Qtt11 / tHsq / tHsqJ
2214  + Quu11 / uHsq / uHsqJ
2215  + Qtu11 / tHsq / uHsqJ;
2216  sigHel[1] += Qtt22 / tHsq / tHsqJ
2217  + Quu22 / uHsq / uHsqJ
2218  + Qtu22 / tHsq / uHsqJ;
2219  sigHel[2] += Qtt12 / tHsq / tHsqJ
2220  + Quu12 / uHsq / uHsqJ
2221  + Qtu12 / tHsq / uHsqJ;
2222  sigHel[3] += Qtt21 / tHsq / tHsqJ
2223  + Quu21 / uHsq / uHsqJ
2224  + Qtu21 / tHsq / uHsqJ;
2225 
2226  }
2227 
2228  }
2229 
2230  // Sum helicity contributions
2231  double sigSum = sigHel[0] + sigHel[1] + sigHel[2] + sigHel[3];
2232 
2233  // Return 0 if all terms vanish, else compute and return cross section
2234  if ( sigSum <= 0. ) return 0.0;
2235 
2236  // Answer
2237  double sigma = (M_PI / 8. / sH2) * pow2(alpS) * sigSum * openFracPair;
2238  return sigma;
2239 
2240 }
2241 
2242 //--------------------------------------------------------------------------
2243 
2244 // Select identity, colour and anticolour.
2245 
2246 void Sigma2qqbar2gluinogluino::setIdColAcol() {
2247 
2248  // Flavours are trivial.
2249  setId( id1, id2, 1000021, 1000021);
2250 
2251  // Two colour flow topologies. Swap if first is antiquark.
2252  if (rndmPtr->flat() < 0.5) setColAcol( 1, 0, 0, 2, 1, 3, 3, 2);
2253  else setColAcol( 1, 0, 0, 2, 3, 2, 1, 3);
2254  if (id1 < 0) swapColAcol();
2255 
2256 }
2257 
2258 //==========================================================================
2259 
2260 // Sigma1qq2antisquark
2261 // R-parity violating squark production
2262 
2263 //--------------------------------------------------------------------------
2264 
2265 // Initialise process
2266 
2267 void Sigma1qq2antisquark::initProc(){
2268 
2269  //Typecast to the correct couplings
2270  coupSUSYPtr = (CoupSUSY*) couplingsPtr;
2271 
2272  //Construct name of the process from lambda'' couplings
2273 
2274  nameSave = "q q' -> " + coupSUSYPtr->getName(idRes)+"* + c.c";
2275  codeSave = 2000 + 10*abs(idRes)/1000000 + abs(idRes)%10;
2276 }
2277 
2278 //--------------------------------------------------------------------------
2279 
2280 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2281 
2282 void Sigma1qq2antisquark::sigmaKin() {
2283 
2284  // Check if at least one RPV coupling non-zero
2285  if(!coupSUSYPtr->isUDD) {
2286  sigBW = 0.0;
2287  return;
2288  }
2289 
2290  mRes = particleDataPtr->m0(abs(idRes));
2291  GammaRes = particleDataPtr->mWidth(abs(idRes));
2292  m2Res = pow2(mRes);
2293 
2294  sigBW = sH * GammaRes/ ( pow2(sH - m2Res) + pow2(mRes * GammaRes) );
2295  sigBW *= 2.0/3.0/mRes;
2296 
2297  // Width out only includes open channels.
2298  widthOut = GammaRes * particleDataPtr->resOpenFrac(id3);
2299 }
2300 
2301 //--------------------------------------------------------------------------
2302 
2303 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2304 
2305 double Sigma1qq2antisquark::sigmaHat() {
2306 
2307  // Only allow (anti)quark-(anti)quark incoming states
2308  if (id1*id2 <= 0) return 0.0;
2309 
2310  // Generation indices
2311  int iA = (abs(id1)+1)/2;
2312  int iB = (abs(id2)+1)/2;
2313 
2314  //Covert from pdg-code to ~u_i/~d_i basis
2315  bool idown = (abs(idRes)%2 == 1) ? true : false;
2316  int iC = (abs(idRes)/1000000 == 2)
2317  ? (abs(idRes)%10+1)/2 + 3: (abs(idRes)%10+1)/2;
2318 
2319  // UDD structure
2320  if (abs(id1)%2 == 0 && abs(id2)%2 == 0) return 0.0;
2321  if (abs(id1)%2 == 1 && abs(id2)%2 == 1 && idown) return 0.0;
2322  if ((abs(id1) + abs(id2))%2 == 1 && !idown) return 0.0;
2323 
2324  double sigma = 0.0;
2325 
2326  if(!idown){
2327  // d_i d_j --> ~u*_k
2328  for(int isq = 1; isq <=3; isq++){
2329  // Loop over R-type squark contributions
2330  sigma += pow2(coupSUSYPtr->rvUDD[isq][iA][iB])
2331  * norm(coupSUSYPtr->Rusq[iC][isq+3]);
2332  }
2333  }else{
2334  // u_i d_j --> d*_k
2335  // Pick the correct coupling for d-u in-state
2336  if(abs(id1)%2==1){
2337  iA = (abs(id2)+1)/2;
2338  iB = (abs(id1)+1)/2;
2339  }
2340  for(int isq = 1; isq <= 3; isq++){
2341  // Loop over R-type squark contributions
2342  sigma += pow2(coupSUSYPtr->rvUDD[iA][iB][isq])
2343  * norm(coupSUSYPtr->Rdsq[iC][isq+3]);
2344  }
2345  }
2346 
2347  sigma *= sigBW;
2348  return sigma;
2349 
2350 }
2351 
2352 //--------------------------------------------------------------------------
2353 
2354 // Select identity, colour and anticolour.
2355 
2356 void Sigma1qq2antisquark::setIdColAcol() {
2357 
2358  // Set flavours.
2359  if(id1 < 0 && id2 < 0 ) setId( id1, id2, idRes);
2360  else setId( id1, id2, -idRes);
2361 
2362  // Colour flow topologies. Swap when antiquarks.
2363  if (abs(id1) < 9) setColAcol( 1, 0, 2, 0, 0, 3);
2364  else setColAcol( 0, 0, 0, 0, 0, 0);
2365  if (id1 < 0) swapColAcol();
2366 
2367 }
2368 
2369 
2370 //==========================================================================
2371 
2372 
2373 // Sigma2qqbar2chi0gluino
2374 // Cross section for gaugino pair production: neutralino-gluino
2375 
2376 //--------------------------------------------------------------------------
2377 
2378 // Initialize process.
2379 
2380 void Sigma2qqbar2chi0gluino::initProc() {
2381 
2382  //Typecast to the correct couplings
2383  coupSUSYPtr = (CoupSUSY*) couplingsPtr;
2384 
2385  // Construct name of process.
2386  nameSave = "q qbar' -> " + particleDataPtr->name(id3) + " "
2387  + particleDataPtr->name(id4);
2388 
2389  // Secondary open width fraction.
2390  openFracPair = particleDataPtr->resOpenFrac(id3, id4);
2391 
2392 }
2393 
2394 //--------------------------------------------------------------------------
2395 
2396 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2397 
2398 void Sigma2qqbar2chi0gluino::sigmaKin() {
2399 
2400  // Common flavour-independent factor.
2401  sigma0 = M_PI * 4.0 / 9.0/ sH2 / coupSUSYPtr->sin2W * alpEM * alpS
2402  * openFracPair;
2403 
2404  // Auxiliary factors for use below
2405  ui = uH - s3;
2406  uj = uH - s4;
2407  ti = tH - s3;
2408  tj = tH - s4;
2409 }
2410 
2411 //--------------------------------------------------------------------------
2412 
2413 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2414 
2415 double Sigma2qqbar2chi0gluino::sigmaHat() {
2416 
2417  // Only allow quark-antiquark incoming states
2418  if (id1*id2 >= 0) return 0.0;
2419 
2420  // In-pair must both be up-type or both down-type
2421  if ((id1+id2) % 2 != 0) return 0.0;
2422 
2423  // Swap T and U if antiquark-quark instead of quark-antiquark
2424  if (id1<0) swapTU = true;
2425 
2426  // Shorthands
2427  int idAbs1 = abs(id1);
2428  int idAbs2 = abs(id2);
2429 
2430  // Flavour-dependent kinematics-dependent couplings.
2431  complex QuLL(0.0),QtLL(0.0),QuRR(0.0),QtRR(0.0);
2432  complex QuLR(0.0),QtLR(0.0),QuRL(0.0),QtRL(0.0);
2433 
2434  // Flavour indices
2435  int ifl1 = (idAbs1+1) / 2;
2436  int ifl2 = (idAbs2+1) / 2;
2437 
2438  complex (*LsddXloc)[4][6];
2439  complex (*RsddXloc)[4][6];
2440  complex (*LsuuXloc)[4][6];
2441  complex (*RsuuXloc)[4][6];
2442  LsddXloc = coupSUSYPtr->LsddX;
2443  RsddXloc = coupSUSYPtr->RsddX;
2444  LsuuXloc = coupSUSYPtr->LsuuX;
2445  RsuuXloc = coupSUSYPtr->RsuuX;
2446 
2447  // Add t-channel squark flavour sums to QmXY couplings
2448  for (int ksq=1; ksq<=6; ksq++) {
2449 
2450  // squark id and squark-subtracted u and t
2451 
2452  int idsq;
2453  idsq=((ksq+2)/3)*1000000 + 2*((ksq-1) % 3) + (idAbs1+1) % 2 + 1;
2454 
2455  double msq2 = pow(particleDataPtr->m0(idsq),2);
2456  double usq = uH - msq2;
2457  double tsq = tH - msq2;
2458 
2459  complex Lsqq1X4;
2460  complex Lsqq2X4;
2461  complex Rsqq1X4;
2462  complex Rsqq2X4;
2463 
2464  complex Lsqq1G;
2465  complex Rsqq1G;
2466  complex Lsqq2G;
2467  complex Rsqq2G;
2468 
2469  // Couplings
2470  Lsqq1X4 = LsuuXloc[ksq][ifl1][id4chi];
2471  Lsqq2X4 = LsuuXloc[ksq][ifl2][id4chi];
2472  Rsqq1X4 = RsuuXloc[ksq][ifl1][id4chi];
2473  Rsqq2X4 = RsuuXloc[ksq][ifl2][id4chi];
2474 
2475  Lsqq1G = coupSUSYPtr->LsuuG[ksq][ifl1];
2476  Rsqq1G = coupSUSYPtr->RsuuG[ksq][ifl1];
2477  Lsqq2G = coupSUSYPtr->LsuuG[ksq][ifl2];
2478  Rsqq2G = coupSUSYPtr->RsuuG[ksq][ifl2];
2479 
2480  if (idAbs1 % 2 != 0) {
2481  Lsqq1X4 = LsddXloc[ksq][ifl1][id4chi];
2482  Lsqq2X4 = LsddXloc[ksq][ifl2][id4chi];
2483  Rsqq1X4 = RsddXloc[ksq][ifl1][id4chi];
2484  Rsqq2X4 = RsddXloc[ksq][ifl2][id4chi];
2485 
2486  Lsqq1G = coupSUSYPtr->LsddG[ksq][ifl1];
2487  Rsqq1G = coupSUSYPtr->RsddG[ksq][ifl1];
2488  Lsqq2G = coupSUSYPtr->LsddG[ksq][ifl2];
2489  Rsqq2G = coupSUSYPtr->RsddG[ksq][ifl2];
2490  }
2491 
2492  // QuXY
2493  QuLL += conj(Lsqq1X4)*Lsqq2G/usq;
2494  QuRR += conj(Rsqq1X4)*Rsqq2G/usq;
2495  QuLR += conj(Lsqq1X4)*Rsqq2G/usq;
2496  QuRL += conj(Rsqq1X4)*Lsqq2G/usq;
2497 
2498  // QtXY
2499  QtLL -= conj(Lsqq1G)*Lsqq2X4/tsq;
2500  QtRR -= conj(Rsqq1G)*Rsqq2X4/tsq;
2501  QtLR += conj(Lsqq1G)*Rsqq2X4/tsq;
2502  QtRL += conj(Rsqq1G)*Lsqq2X4/tsq;
2503 
2504  }
2505 
2506  // Overall factor multiplying coupling
2507  double fac = (1.0-coupSUSYPtr->sin2W);
2508 
2509  // Compute matrix element weight
2510  double weight = 0;
2511  double facLR = uH*tH - s3*s4;
2512  double facMS = m3*m4*sH;
2513 
2514  // Average over separate helicity contributions
2515  // LL (ha = -1, hb = +1) (divided by 4 for average)
2516  weight += norm(QuLL) * ui * uj + norm(QtLL) * ti * tj
2517  + 2 * real(conj(QuLL) * QtLL) * facMS;
2518  // RR (ha = 1, hb = -1) (divided by 4 for average)
2519  weight += norm(QtRR) * ti * tj + norm(QuRR) * ui * uj
2520  + 2 * real(conj(QuRR) * QtRR) * facMS;
2521  // RL (ha = 1, hb = 1) (divided by 4 for average)
2522  weight += norm(QuRL) * ui * uj + norm(QtRL) * ti * tj
2523  + real(conj(QuRL) * QtRL) * facLR;
2524  // LR (ha = -1, hb = -1) (divided by 4 for average)
2525  weight += norm(QuLR) * ui * uj + norm(QtLR) * ti * tj
2526  + real(conj(QuLR) * QtLR) * facLR;
2527 
2528  // Cross section, including colour factor.
2529  double sigma = sigma0 * weight / fac;
2530 
2531  // Answer.
2532  return sigma;
2533 
2534 }
2535 
2536 //--------------------------------------------------------------------------
2537 
2538 // Select identity, colour and anticolour.
2539 
2540 void Sigma2qqbar2chi0gluino::setIdColAcol() {
2541 
2542  // Set flavours.
2543  setId( id1, id2, id3, id4);
2544 
2545  // Colour flow topologies. Swap when antiquarks.
2546  setColAcol( 1, 0, 0, 2, 1, 2, 0, 0);
2547  if (id1 < 0) swapColAcol();
2548 
2549 }
2550 
2551 //==========================================================================
2552 
2553 // Sigma2qqbar2chargluino
2554 // Cross section for gaugino pair production: chargino-gluino
2555 
2556 //--------------------------------------------------------------------------
2557 
2558 // Initialize process.
2559 
2560 void Sigma2qqbar2chargluino::initProc() {
2561 
2562  //Typecast to the correct couplings
2563  coupSUSYPtr = (CoupSUSY*) couplingsPtr;
2564 
2565  // Construct name of process.
2566  nameSave = "q qbar' -> " + particleDataPtr->name(id3) + " "
2567  + particleDataPtr->name(id4) + " + c.c";
2568 
2569  // Secondary open width fraction.
2570  openFracPair = particleDataPtr->resOpenFrac(id3, id4);
2571 
2572 }
2573 
2574 //--------------------------------------------------------------------------
2575 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2576 
2577 void Sigma2qqbar2chargluino::sigmaKin() {
2578 
2579  // Common flavour-independent factor.
2580 
2581  sigma0 = M_PI / sH2 * 4.0 / 9.0 / coupSUSYPtr->sin2W * alpEM * alpS ;
2582  sigma0 /= 2.0 * (1 - coupSUSYPtr->sin2W) ;
2583 
2584  // Auxiliary factors for use below
2585  ui = uH - s3;
2586  uj = uH - s4;
2587  ti = tH - s3;
2588  tj = tH - s4;
2589 }
2590 
2591 //--------------------------------------------------------------------------
2592 
2593 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2594 
2595 double Sigma2qqbar2chargluino::sigmaHat() {
2596 
2597  // Only allow particle-antiparticle incoming states
2598  if (id1*id2 >= 0) return 0.0;
2599 
2600  // Only allow incoming states with sum(charge) = final state
2601  if (abs(id1) % 2 == abs(id2) % 2) return 0.0;
2602  int isPos = (id4chi > 0 ? 1 : 0);
2603  if (id1 < 0 && id1 > -19 && abs(id1) % 2 == 1-isPos ) return 0.0;
2604  else if (id1 > 0 && id1 < 19 && abs(id1) % 2 == isPos ) return 0.0;
2605 
2606  // Flavour-dependent kinematics-dependent couplings.
2607  int idAbs1 = abs(id1);
2608  int iChar = abs(id4chi);
2609 
2610  complex QuLL(0.0),QtLL(0.0),QuRR(0.0),QtRR(0.0);
2611  complex QuLR(0.0),QtLR(0.0),QuRL(0.0),QtRL(0.0);
2612 
2613  // Calculate everything from udbar -> ~chi+ ~chi0 template process
2614  complex LsddGl;
2615  complex RsddGl;
2616  complex LsuuGl;
2617  complex RsuuGl;
2618  complex (*LsduXloc)[4][3];
2619  complex (*RsduXloc)[4][3];
2620  complex (*LsudXloc)[4][3];
2621  complex (*RsudXloc)[4][3];
2622 
2623  LsduXloc = coupSUSYPtr->LsduX;
2624  RsduXloc = coupSUSYPtr->RsduX;
2625  LsudXloc = coupSUSYPtr->LsudX;
2626  RsudXloc = coupSUSYPtr->RsudX;
2627 
2628  // u dbar , ubar d : do nothing
2629  // dbar u , d ubar : swap 1<->2 and t<->u
2630  int iGu = abs(id1)/2;
2631  int iGd = (abs(id2)+1)/2;
2632  if (idAbs1 % 2 != 0) {
2633  swapTU = true;
2634  iGu = abs(id2)/2;
2635  iGd = (abs(id1)+1)/2;
2636  }
2637 
2638  // Add t-channel squark flavour sums to QmXY couplings
2639  for (int jsq=1; jsq<=6; jsq++) {
2640 
2641  int idsu=((jsq+2)/3)*1000000 + 2*((jsq-1) % 3) + 2 ;
2642  int idsd=((jsq+2)/3)*1000000 + 2*((jsq-1) % 3) + 1 ;
2643 
2644  LsddGl = coupSUSYPtr->LsddG[jsq][iGd];
2645  RsddGl = coupSUSYPtr->RsddG[jsq][iGd];
2646  LsuuGl = coupSUSYPtr->LsuuG[jsq][iGu];
2647  RsuuGl = coupSUSYPtr->RsuuG[jsq][iGu];
2648 
2649  double msd2 = pow(particleDataPtr->m0(idsd),2);
2650  double msu2 = pow(particleDataPtr->m0(idsu),2);
2651  double tsq = tH - msd2;
2652  double usq = uH - msu2;
2653 
2654  QuLL += conj(LsuuGl) * conj(LsudXloc[jsq][iGd][iChar])/usq;
2655  QuLR += conj(LsuuGl) * conj(RsudXloc[jsq][iGd][iChar])/usq;
2656  QuRR += conj(RsuuGl) * conj(RsudXloc[jsq][iGd][iChar])/usq;
2657  QuRL += conj(RsuuGl) * conj(LsudXloc[jsq][iGd][iChar])/usq;
2658 
2659  QtLL -= conj(LsduXloc[jsq][iGu][iChar]) * LsddGl/tsq;
2660  QtRR -= conj(RsduXloc[jsq][iGu][iChar]) * RsddGl/tsq;
2661  QtLR += conj(LsduXloc[jsq][iGu][iChar]) * RsddGl/tsq;
2662  QtRL += conj(RsduXloc[jsq][iGu][iChar]) * LsddGl/tsq;
2663  }
2664 
2665  // Compute matrix element weight
2666  double weight = 0;
2667 
2668  // Average over separate helicity contributions
2669  // (if swapped, swap ha, hb if computing polarized cross sections)
2670  // LL (ha = -1, hb = +1) (divided by 4 for average)
2671  weight += norm(QuLL) * ui * uj + norm(QtLL) * ti * tj
2672  + 2 * real(conj(QuLL) * QtLL) * m3 * m4 * sH;
2673  // RR (ha = 1, hb = -1) (divided by 4 for average)
2674  weight += norm(QtRR) * ti * tj + norm(QuRR) * ui * uj
2675  + 2 * real(conj(QuRR) * QtRR) * m3 * m4 * sH;
2676  // RL (ha = 1, hb = 1) (divided by 4 for average)
2677  weight += norm(QuRL) * ui * uj + norm(QtRL) * ti * tj
2678  + real(conj(QuRL) * QtRL) * (uH * tH - s3 * s4);
2679  // LR (ha = -1, hb = -1) (divided by 4 for average)
2680  weight += norm(QuLR) * ui * uj + norm(QtLR) * ti * tj
2681  + real(conj(QuLR) * QtLR) * (uH * tH - s3 * s4);
2682 
2683  // Cross section, including colour factor.
2684  double sigma = sigma0 * weight;
2685 
2686  // Answer.
2687  return sigma;
2688 
2689 }
2690 
2691 //--------------------------------------------------------------------------
2692 
2693 void Sigma2qqbar2chargluino::setIdColAcol() {
2694 
2695  // Set flavours.
2696  setId( id1, id2, id3, id4);
2697 
2698  // Colour flow topologies. Swap when antiquarks.
2699  setColAcol( 1, 0, 0, 2, 1, 2, 0, 0);
2700  if (id1 < 0) swapColAcol();
2701 
2702 }
2703 
2704 //==========================================================================
2705 
2706 // Sigma2qqbar2sleptonantislepton
2707 // Cross section for qqbar-initiated slepton-antislepton production
2708 
2709 //--------------------------------------------------------------------------
2710 
2711 // Initialize process.
2712 
2713 void Sigma2qqbar2sleptonantislepton::initProc() {
2714 
2715  //Typecast to the correct couplings
2716  coupSUSYPtr = (CoupSUSY*) couplingsPtr;
2717 
2718  // Is this a ~u_i ~d*_j, ~d_i ~u*_j final state or ~d_i ~d*_j, ~u_i ~u*_j
2719  if (abs(id3Sav) % 2 == abs(id4Sav) % 2) isUD = false;
2720  else isUD = true;
2721 
2722  // Derive name
2723  nameSave = "q qbar' -> "+particleDataPtr->name(abs(id3Sav))+" "+
2724  particleDataPtr->name(-abs(id4Sav));
2725  if (isUD) nameSave +=" + c.c.";
2726 
2727  // Extract isospin and mass-ordering indices
2728 
2729  if(isUD && abs(id3Sav)%2 == 0) {
2730  // Make sure iGen3 is always slepton and iGen4 is always sneutrino
2731  iGen3 = 3*(abs(id4Sav)/2000000) + (abs(id4Sav)%10+1)/2;
2732  iGen4 = 3*(abs(id3Sav)/2000000) + (abs(id3Sav)%10+1)/2;
2733  }
2734  else {
2735  iGen3 = 3*(abs(id3Sav)/2000000) + (abs(id3Sav)%10+1)/2;
2736  iGen4 = 3*(abs(id4Sav)/2000000) + (abs(id4Sav)%10+1)/2;
2737  }
2738 
2739  // Count 5 neutralinos in NMSSM
2740  nNeut = (coupSUSYPtr->isNMSSM ? 5 : 4);
2741 
2742  // Store mass squares of all possible internal propagator lines;
2743  // retained for future extension to leptonic initial states
2744  m2Neut.resize(nNeut+1);
2745  for (int iNeut=1;iNeut<=nNeut;iNeut++)
2746  m2Neut[iNeut] = pow2(particleDataPtr->m0(coupSUSYPtr->idNeut(iNeut)));
2747 
2748  // Set sizes of some arrays to be used below
2749  tNeut.resize(nNeut+1);
2750  uNeut.resize(nNeut+1);
2751 
2752  // Shorthand for Weak mixing
2753  xW = coupSUSYPtr->sin2W;
2754 
2755  // Secondary open width fraction.
2756  openFracPair = particleDataPtr->resOpenFrac(id3Sav, id4Sav);
2757 
2758 }
2759 
2760 //--------------------------------------------------------------------------
2761 
2762 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
2763 
2764 void Sigma2qqbar2sleptonantislepton::sigmaKin() {
2765 
2766  // Z/W propagator
2767  if (! isUD) {
2768  double sV= sH - pow2(coupSUSYPtr->mZpole);
2769  double d = pow2(sV) + pow2(coupSUSYPtr->mZpole * coupSUSYPtr->wZpole);
2770  propZW = complex( sV / d, coupSUSYPtr->mZpole * coupSUSYPtr->wZpole / d);
2771  } else {
2772  double sV= sH - pow2(coupSUSYPtr->mWpole);
2773  double d = pow2(sV) + pow2(coupSUSYPtr->mWpole * coupSUSYPtr->wWpole);
2774  propZW = complex( sV / d, coupSUSYPtr->mWpole * coupSUSYPtr->wWpole / d);
2775  }
2776 
2777  // Flavor-independent pre-factors
2778  double comFacHat = M_PI/sH2 * openFracPair;
2779 
2780  sigmaEW = comFacHat * pow2(alpEM);
2781 }
2782 
2783 //--------------------------------------------------------------------------
2784 
2785 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
2786 
2787 double Sigma2qqbar2sleptonantislepton::sigmaHat() {
2788 
2789  // In-pair must be opposite-sign
2790  if (id1 * id2 > 0) return 0.0;
2791 
2792  // Check correct charge sum
2793  if (isUD && abs(id1) %2 == abs(id2) % 2) return 0.0;
2794  if (!isUD && abs(id1) % 2 != abs(id2) % 2) return 0.0;
2795 
2796  // No RH sneutrinos
2797  if ( (abs(id3)%2 == 0 && abs(id3) > 2000000)
2798  || (abs(id4)%2 == 0 && abs(id4) > 2000000) ) return 0.0;
2799 
2800  // Coded UD sigma is for udbar -> ~v~l'*. Swap t<->u for dbaru -> ~l~v*.
2801  swapTU = (isUD && abs(id1) % 2 != 0);
2802 
2803  // Coded QQ sigma is for qqbar -> ~l~l*. Swap t<->u for qbarq -> ~l~l*.
2804  if (!isUD && id1 < 0) swapTU = true;
2805 
2806  // Generation indices of incoming particles
2807  int idIn1A = (swapTU) ? abs(id2) : abs(id1);
2808  int idIn2A = (swapTU) ? abs(id1) : abs(id2);
2809  int iGen1 = (idIn1A+1)/2;
2810  int iGen2 = (idIn2A+1)/2;
2811 
2812  // Auxiliary factors for use below
2813  for (int i=1; i<= nNeut; i++) {
2814  tNeut[i] = tH - m2Neut[i];
2815  uNeut[i] = uH - m2Neut[i];
2816  }
2817 
2818  double eQ = (idIn1A % 2 == 0) ? 2./3. : -1./3. ;
2819  double eSl = (abs(id3Sav) % 2 == 0) ? 0. : -1. ;
2820 
2821  // Initial values for pieces used for color-flow selection below
2822  sumColS = 0.0;
2823  sumColT = 0.0;
2824  sumInterference = 0.0;
2825 
2826  // Common factor for LR and RL contributions
2827  double facTU = uH*tH-s3*s4;
2828 
2829  // Opposite-isospin: udbar -> ~l~v*
2830  if ( isUD ) {
2831 
2832  // s-channel W contribution (only contributes to LL helicities)
2833  sumColS = sigmaEW / 32.0 / pow2(xW) / pow2(1.0-xW)
2834  * norm(conj(coupSUSYPtr->LudW[iGen1][iGen2])
2835  * coupSUSYPtr->LslsvW[iGen3][iGen4]) * facTU * norm(propZW);
2836  }
2837 
2838  double CslZ;
2839 
2840  // s-channel Z/photon and interference
2841  if (abs(id1) == abs(id2)) {
2842 
2843  CslZ = real(coupSUSYPtr->LslslZ[iGen3][iGen4]
2844  + coupSUSYPtr->RslslZ[iGen3][iGen4]);
2845  if (abs(id3)%2 == 0)
2846  CslZ = real(coupSUSYPtr->LsvsvZ[iGen3][iGen4]
2847  + coupSUSYPtr->RsvsvZ[iGen3][iGen4]);
2848 
2849  // gamma
2850  // Factor 2 since contributes to both ha != hb helicities
2851  sumColS += (abs(CslZ) > 0.0) ? 2. * pow2(eQ) * pow2(eSl) * sigmaEW
2852  * facTU / pow2(sH) : 0.0;
2853 
2854  // Z/gamma interference
2855  sumInterference += eQ * eSl * sigmaEW * facTU / 2.0 / xW / (1.-xW)
2856  * sqrt(norm(propZW)) / sH * CslZ
2857  * (coupSUSYPtr->LqqZ[idIn1A] + coupSUSYPtr->RqqZ[idIn1A]);
2858 
2859  // s-channel Z
2860 
2861  CslZ = norm(coupSUSYPtr->LslslZ[iGen3][iGen4]
2862  + coupSUSYPtr->RslslZ[iGen3][iGen4]);
2863  if (abs(id3Sav)%2 == 0)
2864  CslZ = norm(coupSUSYPtr->LsvsvZ[iGen3][iGen4]
2865  + coupSUSYPtr->RsvsvZ[iGen3][iGen4]);
2866 
2867  sumColS += sigmaEW * facTU / 16.0 / pow2(xW) / pow2(1.0-xW)
2868  * norm(propZW) * CslZ
2869  * ( pow2(coupSUSYPtr->LqqZ[idIn1A]) + pow2(coupSUSYPtr->RqqZ[idIn1A]) );
2870  }
2871 
2872  // Cross section
2873  double sigma = sumColS + sumColT + sumInterference;
2874 
2875  // Colour average
2876  if(abs(id1) < 10) sigma /= 3.0;
2877 
2878  // Add cc term
2879  if(isUD) sigma *= 2.0;
2880 
2881  // Return answer.
2882  return sigma;
2883 
2884 }
2885 
2886 //--------------------------------------------------------------------------
2887 
2888 // Select identity, colour and anticolour.
2889 
2890 void Sigma2qqbar2sleptonantislepton::setIdColAcol() {
2891 
2892  // Set flavours.
2893  int iSl, iSv;
2894  if( isUD ){
2895  iSl = (abs(id3)%2 == 0) ? abs(id3) : abs(id4);
2896  iSv = (abs(id3)%2 == 0) ? abs(id4) : abs(id3);
2897  if ((id1%2 + id2%2 ) > 0)
2898  setId( id1, id2, -iSl, iSv);
2899  else
2900  setId( id1, id2, iSl, -iSv);
2901  }
2902  else
2903  setId( id1, id2, abs(id3), -abs(id4));
2904 
2905  setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
2906  if (id1 < 0 ) swapColAcol();
2907 
2908 }
2909 
2910 //==========================================================================
2911 
2912 } // end namespace Pythia8
2913 
Event
Definition: AgUStep.h:26
Generated by   doxygen 1.8.5