Update to pythi8.170
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8170 / src / SigmaSUSY.cxx
1 // SigmaSUSY.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2012 Torbjorn Sjostrand.
3 // Main authors of this file: N. Desai, P. Skands
4 // PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
5 // Please respect the MCnet Guidelines, see GUIDELINES for details.
6
7 // Function definitions (not found in the header) for the 
8 // supersymmetry simulation classes. 
9
10 #include "SigmaSUSY.h"
11
12 namespace Pythia8 {
13   
14 //==========================================================================
15
16 // Sigma2qqbar2chi0chi0 
17 // Cross section for gaugino pair production: neutralino pair
18
19 //--------------------------------------------------------------------------
20
21 // Initialize process. 
22   
23 void Sigma2qqbar2chi0chi0::initProc() {
24
25   //Typecast to the correct couplings
26   coupSUSYPtr = (CoupSUSY*) couplingsPtr;
27
28   // Construct name of process. 
29   nameSave = "q qbar' -> " + particleDataPtr->name(id3) + " " 
30     + particleDataPtr->name(id4);
31
32   // Secondary open width fraction.
33   openFracPair = particleDataPtr->resOpenFrac(id3, id4);
34
35 }
36
37 //--------------------------------------------------------------------------
38
39 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour. 
40
41 void Sigma2qqbar2chi0chi0::sigmaKin() {
42
43   // Common flavour-independent factor.
44   sigma0 = M_PI /3.0/ sH2 / pow2(coupSUSYPtr->sin2W) * pow2(alpEM) 
45     * openFracPair; 
46
47   // Auxiliary factors for use below
48   ui       = uH - s3;
49   uj       = uH - s4;
50   ti       = tH - s3;
51   tj       = tH - s4;
52   double sV= sH - pow2(coupSUSYPtr->mZpole);
53   double d = pow2(sV) + pow2(coupSUSYPtr->mZpole * coupSUSYPtr->wZpole);
54   propZ    = complex( sV / d, coupSUSYPtr->mZpole * coupSUSYPtr->wZpole / d);
55
56 }
57
58 //--------------------------------------------------------------------------
59
60 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence. 
61
62 double Sigma2qqbar2chi0chi0::sigmaHat() {
63
64   // Only allow quark-antiquark incoming states
65   if (id1*id2 >= 0) {
66     return 0.0;    
67   }
68   
69   // Only allow incoming states with sum(charge) = 0
70   if ((id1+id2) % 2 != 0) {
71     return 0.0;    
72   }
73
74   if(id1<0) swapTU = true;
75   
76   // Shorthands
77   int idAbs1    = abs(id1);  
78   int idAbs2    = abs(id2);
79
80   // Flavour-dependent kinematics-dependent couplings.
81   complex QuLL(0.0),QtLL(0.0),QuRR(0.0),QtRR(0.0);
82   complex QuLR(0.0),QtLR(0.0),QuRL(0.0),QtRL(0.0);
83
84   // s-channel Z couplings
85   if (idAbs1 == idAbs2) {
86     QuLL = coupSUSYPtr->LqqZ[idAbs1] * coupSUSYPtr->OLpp[id3chi][id4chi] 
87          * propZ / 2.0;
88     QtLL = coupSUSYPtr->LqqZ[idAbs1] * coupSUSYPtr->ORpp[id3chi][id4chi] 
89          * propZ / 2.0;
90     QuRR = coupSUSYPtr->RqqZ[idAbs1] * coupSUSYPtr->ORpp[id3chi][id4chi] 
91          * propZ / 2.0;
92     QtRR = coupSUSYPtr->RqqZ[idAbs1] * coupSUSYPtr->OLpp[id3chi][id4chi] 
93          * propZ / 2.0;
94   }
95
96   // Flavour indices
97   int ifl1 = (idAbs1+1) / 2;
98   int ifl2 = (idAbs2+1) / 2;
99
100   // Add t-channel squark flavour sums to QmXY couplings
101   for (int ksq=1; ksq<=6; ksq++) {    
102
103     // squark id and squark-subtracted u and t
104     int idsq=((ksq+2)/3)*1000000 + 2*((ksq-1) % 3) + (idAbs1+1) % 2 + 1;
105     double msq2    = pow(particleDataPtr->m0(idsq),2);
106     double usq     = uH - msq2;
107     double tsq     = tH - msq2;
108     
109     // Couplings
110     complex Lsqq1X3 = coupSUSYPtr->LsuuX[ksq][ifl1][id3chi];
111     complex Lsqq1X4 = coupSUSYPtr->LsuuX[ksq][ifl1][id4chi];
112     complex Lsqq2X3 = coupSUSYPtr->LsuuX[ksq][ifl2][id3chi];
113     complex Lsqq2X4 = coupSUSYPtr->LsuuX[ksq][ifl2][id4chi];
114     complex Rsqq1X3 = coupSUSYPtr->RsuuX[ksq][ifl1][id3chi];
115     complex Rsqq1X4 = coupSUSYPtr->RsuuX[ksq][ifl1][id4chi];
116     complex Rsqq2X3 = coupSUSYPtr->RsuuX[ksq][ifl2][id3chi];
117     complex Rsqq2X4 = coupSUSYPtr->RsuuX[ksq][ifl2][id4chi];
118     if (idAbs1 % 2 != 0) {
119       Lsqq1X3 = coupSUSYPtr->LsddX[ksq][ifl1][id3chi];
120       Lsqq1X4 = coupSUSYPtr->LsddX[ksq][ifl1][id4chi];
121       Lsqq2X3 = coupSUSYPtr->LsddX[ksq][ifl2][id3chi];
122       Lsqq2X4 = coupSUSYPtr->LsddX[ksq][ifl2][id4chi];
123       Rsqq1X3 = coupSUSYPtr->RsddX[ksq][ifl1][id3chi];
124       Rsqq1X4 = coupSUSYPtr->RsddX[ksq][ifl1][id4chi];
125       Rsqq2X3 = coupSUSYPtr->RsddX[ksq][ifl2][id3chi];
126       Rsqq2X4 = coupSUSYPtr->RsddX[ksq][ifl2][id4chi];      
127     }
128
129     // QuXY
130     QuLL += conj(Lsqq1X4)*Lsqq2X3/usq;
131     QuRR += conj(Rsqq1X4)*Rsqq2X3/usq;
132     QuLR += conj(Lsqq1X4)*Rsqq2X3/usq;
133     QuRL += conj(Rsqq1X4)*Lsqq2X3/usq;
134     
135     
136     // QtXY
137     QtLL -= conj(Lsqq1X3)*Lsqq2X4/tsq;
138     QtRR -= conj(Rsqq1X3)*Rsqq2X4/tsq;
139     QtLR += conj(Lsqq1X3)*Rsqq2X4/tsq;
140     QtRL += conj(Rsqq1X3)*Lsqq2X4/tsq;
141
142   }
143
144   // Overall factor multiplying each coupling; multiplied at the end as fac^2
145   double fac = (1.0-coupSUSYPtr->sin2W);
146   if(abs(id3)==abs(id4)) fac *= sqrt(2.); // for identical final particles
147
148   // Compute matrix element weight
149   double weight = 0;
150   double facLR = uH*tH - s3*s4;
151   double facMS = m3*m4*sH;
152
153   // Average over separate helicity contributions
154   // LL (ha = -1, hb = +1) (divided by 4 for average)            
155   weight += norm(QuLL) * ui * uj + norm(QtLL) * ti * tj
156     + 2 * real(conj(QuLL) * QtLL) * facMS;
157   // RR (ha =  1, hb = -1) (divided by 4 for average)        
158   weight += norm(QtRR) * ti * tj + norm(QuRR) * ui * uj  
159     + 2 * real(conj(QuRR) * QtRR) * facMS;
160   // RL (ha =  1, hb =  1) (divided by 4 for average)        
161   weight += norm(QuRL) * ui * uj + norm(QtRL) * ti * tj
162     + real(conj(QuRL) * QtRL) * facLR;
163   // LR (ha = -1, hb = -1) (divided by 4 for average)        
164   weight += norm(QuLR) * ui * uj + norm(QtLR) * ti * tj
165     + real(conj(QuLR) * QtLR) * facLR;
166
167   // Cross section, including colour factor.
168   double sigma = sigma0 * weight / pow2(fac);
169
170   // Answer.
171   return sigma;    
172
173 }
174
175 //--------------------------------------------------------------------------
176
177 // Select identity, colour and anticolour.
178
179 void Sigma2qqbar2chi0chi0::setIdColAcol() {
180
181   // Set flavours.
182   setId( id1, id2, id3, id4);
183
184   // Colour flow topologies. Swap when antiquarks.
185   if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
186   else              setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
187   if (id1 < 0) swapColAcol();
188
189 }
190
191 //==========================================================================
192
193 // Sigma2qqbar2charchi0
194 // Cross section for gaugino pair production: neutralino-chargino
195
196 //--------------------------------------------------------------------------
197
198 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour. 
199
200 void Sigma2qqbar2charchi0::sigmaKin() {
201
202   // Common flavour-independent factor.
203   
204   sigma0 = M_PI / sH2 / 3.0 / pow2(coupSUSYPtr->sin2W) * pow2(alpEM) ; 
205   sigma0 /= 2.0 * (1 - coupSUSYPtr->sin2W) ; 
206
207   // Auxiliary factors for use below
208   ui        = uH - s3;
209   uj        = uH - s4;
210   ti        = tH - s3;
211   tj        = tH - s4;
212   double sW = sH - pow2(coupSUSYPtr->mWpole);
213   double d  = pow2(sW) + pow2(coupSUSYPtr->mWpole * coupSUSYPtr->wWpole);
214   propW     = complex( sW / d, coupSUSYPtr->mWpole * coupSUSYPtr->wWpole / d);
215
216 }
217
218 //--------------------------------------------------------------------------
219
220 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence. 
221
222 double Sigma2qqbar2charchi0::sigmaHat() {
223
224   // Only allow particle-antiparticle incoming states
225   if (id1*id2 >= 0) {
226     return 0.0;    
227   }
228   
229   // Only allow incoming states with sum(charge) = final state
230   if (abs(id1) % 2 == abs(id2) % 2) return 0.0;
231   int isPos  = (id3chi > 0 ? 1 : 0);
232   if (id1 < 0 && id1 > -10 && abs(id1) % 2 == 1-isPos ) return 0.0;
233   else if (id1 > 0 && id1 < 10 && abs(id1) % 2 == isPos ) return 0.0;
234
235   // Flavour-dependent kinematics-dependent couplings.
236   int idAbs1  = abs(id1);  
237   int iChar = abs(id3chi);
238   int iNeut = abs(id4chi); 
239   
240   complex QuLL(0.0),QtLL(0.0),QuRR(0.0),QtRR(0.0);
241   complex QuLR(0.0),QtLR(0.0),QuRL(0.0),QtRL(0.0);
242   
243   // Calculate everything from udbar -> ~chi+ ~chi0 template process
244   
245   // u dbar , ubar d : do nothing
246   // dbar u , d ubar : swap 1<->2 and t<->u
247
248   int iGu = abs(id1)/2;
249   int iGd = (abs(id2)+1)/2;
250
251   if (idAbs1 % 2 != 0) {
252     swapTU = true;
253     iGu = abs(id2)/2;
254     iGd = (abs(id1)+1)/2;
255   }
256
257   // s-channel W contribution
258   QuLL = conj(coupSUSYPtr->LudW[iGu][iGd]) 
259     * conj(coupSUSYPtr->OL[iNeut][iChar])
260     * propW / sqrt(2.0);
261   QtLL = conj(coupSUSYPtr->LudW[iGu][iGd]) 
262     * conj(coupSUSYPtr->OR[iNeut][iChar])
263     * propW / sqrt(2.0);
264
265   // Add t-channel squark flavour sums to QmXY couplings
266   for (int jsq=1; jsq<=6; jsq++) {    
267
268     int idsu=((jsq+2)/3)*1000000 + 2*((jsq-1) % 3) + 2;
269     int idsd=((jsq+2)/3)*1000000 + 2*((jsq-1) % 3) + 1;
270     double msd2 = pow(particleDataPtr->m0(idsd),2);
271     double msu2 = pow(particleDataPtr->m0(idsu),2);
272     double tsq  = tH - msd2;
273     double usq  = uH - msu2;
274
275     QuLL += conj(coupSUSYPtr->LsuuX[jsq][iGu][iNeut])
276       * conj(coupSUSYPtr->LsudX[jsq][iGd][iChar])/usq;
277     QuLR += conj(coupSUSYPtr->LsuuX[jsq][iGu][iNeut])
278       * conj(coupSUSYPtr->RsudX[jsq][iGd][iChar])/usq;
279     QuRR += conj(coupSUSYPtr->RsuuX[jsq][iGu][iNeut])
280       * conj(coupSUSYPtr->RsudX[jsq][iGd][iChar])/usq;
281     QuRL += conj(coupSUSYPtr->RsuuX[jsq][iGu][iNeut])
282       * conj(coupSUSYPtr->LsudX[jsq][iGd][iChar])/usq;
283
284     QtLL -= conj(coupSUSYPtr->LsduX[jsq][iGu][iChar])
285       * coupSUSYPtr->LsddX[jsq][iGd][iNeut]/tsq;
286     QtRR -= conj(coupSUSYPtr->RsduX[jsq][iGu][iChar])
287       * coupSUSYPtr->RsddX[jsq][iGd][iNeut]/tsq;
288     QtLR += conj(coupSUSYPtr->LsduX[jsq][iGu][iChar])
289       * coupSUSYPtr->RsddX[jsq][iGd][iNeut]/tsq;
290     QtRL += conj(coupSUSYPtr->RsduX[jsq][iGu][iChar])
291       * coupSUSYPtr->LsddX[jsq][iGd][iNeut]/tsq;
292   }
293
294   // Compute matrix element weight
295   double weight = 0;
296
297   // Average over separate helicity contributions
298   // (if swapped, swap ha, hb if computing polarized cross sections)
299   // LL (ha = -1, hb = +1) (divided by 4 for average)            
300   weight += norm(QuLL) * ui * uj + norm(QtLL) * ti * tj
301     + 2 * real(conj(QuLL) * QtLL) * m3 * m4 * sH;
302   // RR (ha =  1, hb = -1) (divided by 4 for average)        
303   weight += norm(QtRR) * ti * tj + norm(QuRR) * ui * uj  
304     + 2 * real(conj(QuRR) * QtRR) * m3 * m4 * sH;
305   // RL (ha =  1, hb =  1) (divided by 4 for average)        
306   weight += norm(QuRL) * ui * uj + norm(QtRL) * ti * tj
307     + real(conj(QuRL) * QtRL) * (uH * tH - s3 * s4);
308   // LR (ha = -1, hb = -1) (divided by 4 for average)        
309   weight += norm(QuLR) * ui * uj + norm(QtLR) * ti * tj
310     + real(conj(QuLR) * QtLR) * (uH * tH - s3 * s4);
311
312   // Cross section, including colour factor.
313   double sigma = sigma0 * weight;
314
315   // Answer.
316   return sigma;    
317
318 }
319
320 //==========================================================================
321
322 // Sigma2qqbar2charchar
323 // Cross section for gaugino pair production: chargino-chargino
324
325 //--------------------------------------------------------------------------
326
327 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour. 
328
329 void Sigma2qqbar2charchar::sigmaKin() {
330
331   // Common flavour-independent factor.
332   sigma0 = M_PI / 3.0 / sH2 / pow2(coupSUSYPtr->sin2W) * pow2(alpEM) ; 
333
334   // Auxiliary factors for use below
335   ui       = uH - s3;
336   uj       = uH - s4;
337   ti       = tH - s3;
338   tj       = tH - s4;
339   double sV= sH - pow2(coupSUSYPtr->mZpole);
340   double d = pow2(sV) + pow2(coupSUSYPtr->mZpole * coupSUSYPtr->wZpole);
341   propZ    = complex( sV / d, coupSUSYPtr->mZpole * coupSUSYPtr->wZpole / d);
342
343 }
344
345 //--------------------------------------------------------------------------
346
347 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence. 
348
349 double Sigma2qqbar2charchar::sigmaHat() { 
350
351   // Only allow quark-antiquark incoming states
352   if (id1*id2 >= 0) {
353     return 0.0;
354   }
355   
356   // Only allow incoming states with sum(charge) = 0
357   if ((id1+id2) % 2 != 0) {
358     return 0.0;    
359   }
360   
361   //if (id1 > 0 || id1==-1 || id1==-3 || id1==-5) return 0.0;
362   //if (id1 < 0 || id1==1 || id1==3 || id1==5) return 0.0;
363   
364   swapTU = (id1 < 0 ? true : false);
365
366   // Flavour-dependent kinematics-dependent couplings.
367   int idAbs1    = abs(id1);  
368   int idAbs2    = abs(id2);  
369   int i3        = abs(id3chi);
370   int i4        = abs(id4chi);
371   
372   // Flavour-dependent kinematics-dependent couplings.
373   complex QuLL(0.0),QtLL(0.0),QuRR(0.0),QtRR(0.0);
374   complex QuLR(0.0),QtLR(0.0),QuRL(0.0),QtRL(0.0);
375
376   // Add Z/gamma* for same-flavour in-quarks
377   if (idAbs1 == idAbs2) {
378     
379     QuLL = -coupSUSYPtr->LqqZ[idAbs1]*conj(coupSUSYPtr->ORp[i3][i4]);
380     QtLL = -coupSUSYPtr->LqqZ[idAbs1]*conj(coupSUSYPtr->OLp[i3][i4]);
381     QuRR = -coupSUSYPtr->RqqZ[idAbs1]*conj(coupSUSYPtr->OLp[i3][i4]);
382     QtRR = -coupSUSYPtr->RqqZ[idAbs1]*conj(coupSUSYPtr->ORp[i3][i4]);
383
384     QuLL *= propZ / 2.0 / (1.0-coupSUSYPtr->sin2W);
385     QtLL *= propZ / 2.0 / (1.0-coupSUSYPtr->sin2W);
386     QuRR *= propZ / 2.0 / (1.0-coupSUSYPtr->sin2W);
387     QtRR *= propZ / 2.0 / (1.0-coupSUSYPtr->sin2W);  
388   
389     // s-channel gamma* (only for same-type charginos)
390     if (i3 == i4) {
391     
392       // Charge of in-particles
393       double q = 2.0/3.0;
394       if (idAbs1 % 2 == 1) q = -1.0/3.0;      
395       QuLL += q * coupSUSYPtr->sin2W / sH;
396       QuRR += q * coupSUSYPtr->sin2W / sH;
397       QtLL += q * coupSUSYPtr->sin2W / sH;
398       QtRR += q * coupSUSYPtr->sin2W / sH;
399     
400     }
401   }
402
403   int iG1    = (abs(id1)+1)/2;
404   int iG2    = (abs(id2)+1)/2;
405
406   // Add t- or u-channel squark flavour sums to QmXY couplings
407   for (int ksq=1; ksq<=6; ksq++) {    
408     
409     if(id1 % 2 == 0) { 
410
411       // u-channel diagrams only
412       // up-type incoming; u-channel ~d 
413       
414       int idsd    = ((ksq+2)/3)*1000000 + 2*((ksq-1) % 3) + 1;
415       double msq  = particleDataPtr->m0(idsd);
416       double ufac = 2.0 * (uH - pow2(msq));
417
418       //u-ubar -> chi-chi+
419       QuLL += coupSUSYPtr->LsduX[ksq][iG2][i3] 
420             * conj(coupSUSYPtr->LsduX[ksq][iG1][i4]) / ufac;
421       QuRR += coupSUSYPtr->RsduX[ksq][iG2][i3] 
422             * conj(coupSUSYPtr->RsduX[ksq][iG1][i4]) / ufac;
423       QuLR += coupSUSYPtr->RsduX[ksq][iG2][i3] 
424             * conj(coupSUSYPtr->LsduX[ksq][iG1][i4]) / ufac;
425       QuRL += coupSUSYPtr->LsduX[ksq][iG2][i3] 
426             * conj(coupSUSYPtr->RsduX[ksq][iG1][i4]) / ufac;
427
428
429     }else{
430       // t-channel diagrams only;
431       // down-type incoming; t-channel ~u 
432
433       int idsu    = ((ksq+2)/3)*1000000 + 2*((ksq-1) % 3) + 2;
434       double msq  = particleDataPtr->m0(idsu);
435       double tfac = 2.0 * (tH - pow2(msq));
436
437       //d-dbar -> chi-chi+
438       QtLL -= coupSUSYPtr->LsudX[ksq][iG1][i3] 
439             * conj(coupSUSYPtr->LsudX[ksq][iG2][i4]) / tfac;
440       QtRR -= coupSUSYPtr->RsudX[ksq][iG1][i3] 
441             * conj(coupSUSYPtr->RsudX[ksq][iG2][i4]) / tfac;
442       QtLR += coupSUSYPtr->LsudX[ksq][iG1][i3] 
443             * conj(coupSUSYPtr->RsudX[ksq][iG2][i4]) / tfac;
444       QtRL += coupSUSYPtr->RsudX[ksq][iG1][i3] 
445             * conj(coupSUSYPtr->LsudX[ksq][iG2][i4]) / tfac;
446
447     }
448   }
449    // Compute matrix element weight
450    double weight = 0;
451
452   // Average over separate helicity contributions
453   // LL (ha = -1, hb = +1) (divided by 4 for average)            
454   weight += norm(QuLL) * ui * uj + norm(QtLL) * ti * tj
455     + 2 * real(conj(QuLL) * QtLL) * m3 * m4 * sH;
456   // RR (ha =  1, hb = -1) (divided by 4 for average)        
457   weight += norm(QtRR) * ti * tj + norm(QuRR) * ui * uj  
458     + 2 * real(conj(QuRR) * QtRR) * m3 * m4 * sH;
459   // RL (ha =  1, hb =  1) (divided by 4 for average)        
460   weight += norm(QuRL) * ui * uj + norm(QtRL) * ti * tj
461     + real(conj(QuRL) * QtRL) * (uH * tH - s3 * s4);
462   // LR (ha = -1, hb = -1) (divided by 4 for average)        
463   weight += norm(QuLR) * ui * uj + norm(QtLR) * ti * tj
464     + real(conj(QuLR) * QtLR) * (uH * tH - s3 * s4);
465
466   // Cross section, including colour factor.
467   double sigma = sigma0 * weight;
468
469   // Answer.
470   return sigma;    
471
472 }
473
474 //==========================================================================
475
476 // Sigma2qgchi0squark 
477 // Cross section for gaugino-squark production: neutralino-squark
478
479 //--------------------------------------------------------------------------
480
481 // Initialize process. 
482   
483 void Sigma2qg2chi0squark::initProc() {
484
485   //Typecast to the correct couplings
486   coupSUSYPtr = (CoupSUSY*) couplingsPtr;
487
488   // Construct name of process. 
489   if (id4 % 2 == 0) {
490     nameSave = "q g -> " + particleDataPtr->name(id3) + " " 
491       + particleDataPtr->name(id4) + " + c.c. (q=u,c)";
492   } 
493   else {
494     nameSave = "q g -> " + particleDataPtr->name(id3) + " " 
495       + particleDataPtr->name(id4) + " + c.c. (q=d,s,b)";
496   }
497
498   // Secondary open width fraction.
499   openFracPair = particleDataPtr->resOpenFrac(id3, id4);
500
501 }
502
503 //--------------------------------------------------------------------------
504
505 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour. 
506
507 void Sigma2qg2chi0squark::sigmaKin() {
508
509   // Common flavour-independent factor.
510   // tmp: alphaS = 0.1 for counter-checks
511   sigma0 = M_PI / sH2 / coupSUSYPtr->sin2W * alpEM * alpS
512     * openFracPair;
513
514   // Auxiliary factors for use below
515   ui       = uH - s3;
516   uj       = uH - s4;
517   ti       = tH - s3;
518   tj       = tH - s4;
519
520 }
521
522 //--------------------------------------------------------------------------
523
524 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence. 
525
526 double Sigma2qg2chi0squark::sigmaHat() {
527
528   // Antiquark -> antisquark
529   int idq = id1;
530   if (id1 == 21 || id1 == 22) idq = id2;
531   if (idq < 0) {
532     id4 = -abs(id4);
533   } else {
534     id4 = abs(id4);
535   }
536
537   // tmp: only allow incoming quarks on side 1
538   //  if (id1 < 0 || id1 == 21) return 0.0;
539
540   // Generation index
541   int iGq = (abs(idq)+1)/2;
542
543   // Only accept u(bar) -> ~u(bar) and d(bar) -> ~d(bar)
544   if (particleDataPtr->chargeType(idq) != particleDataPtr->chargeType(id4))
545     return 0.0;
546   
547   // Couplings
548   complex LsqqX, RsqqX;
549   if (idq % 2 == 0) {
550     LsqqX = coupSUSYPtr->LsuuX[id4sq][iGq][id3chi];
551     RsqqX = coupSUSYPtr->RsuuX[id4sq][iGq][id3chi];
552   }
553   else { 
554     LsqqX = coupSUSYPtr->LsddX[id4sq][iGq][id3chi];
555     RsqqX = coupSUSYPtr->RsddX[id4sq][iGq][id3chi];
556   }  
557
558   // Prefactors : swap u and t if gq instead of qg
559   double fac1, fac2;
560   if (idq == id1) {
561     fac1 = -ui/sH + 2.0 * ( uH*tH - s4*s3 )/sH/tj;
562     fac2 = ti/tj * ( (tH + s4)/tj + (ti - uj)/sH );
563   } else {
564     fac1 = -ti/sH + 2.0 * ( uH*tH - s4*s3 )/sH/uj;
565     fac2 = ui/uj * ( (uH + s4)/uj + (ui - tj)/sH );
566   }
567
568   // Compute matrix element weight
569   double weight = 0.0;
570
571   // Average over separate helicity contributions
572   // (for qbar g : ha -> -ha )
573   // LL (ha = -1, hb = +1) (divided by 4 for average)            
574   weight += fac2 * norm(LsqqX) / 2.0;
575   // RR (ha =  1, hb = -1) (divided by 4 for average)        
576   weight += fac2 * norm(RsqqX) / 2.0;
577   // RL (ha =  1, hb =  1) (divided by 4 for average)        
578   weight += fac2 * norm(RsqqX) / 2.0 + fac1 * norm(RsqqX);
579   // LR (ha = -1, hb = -1) (divided by 4 for average)        
580   weight += fac2 * norm(LsqqX) / 2.0 + fac1 * norm(LsqqX);
581
582   double sigma = sigma0 * weight;
583   if (abs(idq) < 9) sigma /= 3.;
584
585   // Answer.
586   return sigma;    
587
588 }
589
590 //--------------------------------------------------------------------------
591
592 // Select identity, colour and anticolour.
593
594 void Sigma2qg2chi0squark::setIdColAcol() {
595
596   // Set flavours.
597   setId( id1, id2, id3, (id1*id2 > 0 ? abs(id4) : -abs(id4)));
598
599   // Colour flow topology. Swap if first is gluon, or when antiquark.
600   if (id1 != 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
601   else setColAcol( 1, 2, 2, 0, 0, 0, 1, 0);
602   if (id1*id2 < 0) swapColAcol();
603
604 }
605
606 //==========================================================================
607
608 // Sigma2qg2charsquark
609 // Cross section for gaugino-squark production: chargino-squark
610
611 //--------------------------------------------------------------------------
612
613 // Initialize process. 
614   
615 void Sigma2qg2charsquark::initProc() {
616
617   //Typecast to the correct couplings
618   coupSUSYPtr = (CoupSUSY*) couplingsPtr;
619
620   // Construct name of process. 
621   if (id4 % 2 == 0) {
622     nameSave = "q g -> " + particleDataPtr->name(id3) + " " 
623       + particleDataPtr->name(id4) + " + c.c. (q=d,s,b)";
624   } 
625   else {
626     nameSave = "q g -> " + particleDataPtr->name(id3) + " " 
627       + particleDataPtr->name(id4) + " + c.c. (q=u,c)";
628   }
629
630   // Secondary open width fraction.
631   openFracPair = particleDataPtr->resOpenFrac(id3Sav, id4Sav);
632
633 }
634
635 //--------------------------------------------------------------------------
636
637 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence. 
638
639 double Sigma2qg2charsquark::sigmaHat() {
640
641   // Antiquark -> antisquark
642   int idq = (id1 == 21) ? id2 : id1;
643   if (idq > 0) {
644     id3 = id3Sav;
645     id4 = id4Sav;
646   } else {
647     id3 = -id3Sav;
648     id4 = -id4Sav;
649   }
650
651   // Only accept u(bar) -> ~d(bar) and d(bar) -> ~u(bar)
652   if (particleDataPtr->chargeType(idq) == particleDataPtr->chargeType(id4))
653     return 0.0;
654   
655   // Generation index
656   int iGq = (abs(idq)+1)/2;
657
658   // Couplings
659   complex LsqqX, RsqqX;
660   if (idq % 2 == 0) {
661     LsqqX = coupSUSYPtr->LsduX[id4sq][iGq][id3chi];
662     RsqqX = coupSUSYPtr->RsduX[id4sq][iGq][id3chi];
663   }
664   else { 
665     LsqqX = coupSUSYPtr->LsduX[id4sq][iGq][id3chi];
666     RsqqX = coupSUSYPtr->RsduX[id4sq][iGq][id3chi];
667   }  
668
669   // Prefactors : swap u and t if gq instead of qg
670   double fac1, fac2;
671   if (idq == id1) {
672     fac1 = -ui/sH + 2.0 * ( uH*tH - s4*s3 )/sH/tj;
673     fac2 = ti/tj * ( (tH + s4)/tj + (ti - uj)/sH );
674   } else {
675     fac1 = -ti/sH + 2.0 * ( uH*tH - s4*s3 )/sH/uj;
676     fac2 = ui/uj * ( (uH + s4)/uj + (ui - tj)/sH );
677   }
678
679   // Compute matrix element weight
680   double weight = 0.0;
681
682   // Average over separate helicity contributions
683   // (a, b refers to qg configuration)
684   // LL (ha = -1, hb = +1) (divided by 4 for average)            
685   weight += fac2 * norm(LsqqX) / 2.0;
686   // RR (ha =  1, hb = -1) (divided by 4 for average)        
687   weight += fac2 * norm(RsqqX) / 2.0;
688   // RL (ha =  1, hb =  1) (divided by 4 for average)        
689   weight += fac2 * norm(RsqqX) / 2.0 + fac1 * norm(RsqqX);
690   // LR (ha = -1, hb = -1) (divided by 4 for average)        
691   weight += fac2 * norm(LsqqX) / 2.0 + fac1 * norm(LsqqX);
692
693   double sigma = sigma0 * weight;
694   if (abs(idq) < 9) sigma /= 3.;
695
696   // Answer.
697   return sigma * openFracPair;    
698
699 }
700
701 //--------------------------------------------------------------------------
702
703 // Select identity, colour and anticolour.
704
705 void Sigma2qg2charsquark::setIdColAcol() {
706
707   // Set flavours.
708   if (id1 > 0 && id2 > 0) {
709     setId( id1, id2, id3Sav, id4Sav);
710   } else {
711     setId( id1, id2,-id3Sav,-id4Sav);
712   }
713
714   // Colour flow topology. Swap if first is gluon, or when antiquark.
715   if (id1 != 21) setColAcol( 1, 0, 2, 1, 0, 0, 2, 0);
716   else setColAcol( 1, 2, 2, 0, 0, 0, 1, 0);
717   if (id1 < 0 || id2 < 0) swapColAcol();
718
719 }
720
721 //==========================================================================
722
723 // Sigma2qq2squarksquark
724 // Cross section for squark-squark production
725
726 //--------------------------------------------------------------------------
727
728 // Initialize process. 
729   
730 void Sigma2qq2squarksquark::initProc() {
731
732   //Typecast to the correct couplings
733   coupSUSYPtr = (CoupSUSY*) couplingsPtr;
734
735   // Extract mass-ordering indices
736   iGen3 = 3*(abs(id3Sav)/2000000) + (abs(id3Sav)%10+1)/2;
737   iGen4 = 3*(abs(id4Sav)/2000000) + (abs(id4Sav)%10+1)/2;
738
739   // Is this a ~u_i ~d_j fial state or ~d_i ~d_j, ~u_i ~u_j
740   if (abs(id3Sav) % 2 == abs(id4Sav) % 2) isUD = false;
741   else isUD = true;
742
743   // Derive name
744   nameSave = "q q' -> "+particleDataPtr->name(abs(id3Sav))+" "
745     +particleDataPtr->name(abs(id4Sav))+" + c.c.";
746
747   // Count 5 neutralinos in NMSSM
748   nNeut = (coupSUSYPtr->isNMSSM ? 5 : 4);
749
750   // Store mass squares of all possible internal propagator lines
751   m2Glu     = pow2(particleDataPtr->m0(1000021));
752   m2Neut.resize(nNeut+1);
753   for (int iNeut=1;iNeut<=nNeut;iNeut++) {
754     m2Neut[iNeut] = pow2(particleDataPtr->m0(coupSUSYPtr->idNeut(iNeut)));
755   }
756   m2Char.resize(3);
757   m2Char[1] = pow2(particleDataPtr->m0(coupSUSYPtr->idChar(1)));
758   m2Char[2] = pow2(particleDataPtr->m0(coupSUSYPtr->idChar(2)));
759   
760   // Set sizes of some arrays to be used below
761   tNeut.resize(nNeut+1);
762   uNeut.resize(nNeut+1);
763   tChar.resize(3);
764   uChar.resize(3);
765
766   // Secondary open width fraction.
767   openFracPair = particleDataPtr->resOpenFrac(id3Sav, id4Sav);
768 }
769
770 //--------------------------------------------------------------------------
771
772 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour. 
773
774 void Sigma2qq2squarksquark::sigmaKin() {
775
776   // Weak mixing
777   double xW = coupSUSYPtr->sin2W;
778
779   // pi/sH2
780   double comFacHat = M_PI/sH2 * openFracPair;
781
782   // Channel-dependent but flavor-independent pre-factors
783   sigmaNeut     = comFacHat * pow2(alpEM) / pow2(xW) / pow2(1-xW);
784   sigmaGlu      = comFacHat * 2.0 * pow2(alpS) / 9.0;
785   if (isUD) {
786     sigmaChar     = comFacHat * pow2(alpEM) / 4.0 / pow2(xW);
787     sigmaCharNeut = comFacHat * pow2(alpEM) / 3.0 / pow2(xW) / (1-xW); 
788     sigmaCharGlu  = comFacHat * 4.0 * alpEM * alpS / 9.0 / xW;
789     sigmaNeutGlu  = 0.0;
790   } else {
791     sigmaChar     = 0.0;
792     sigmaCharNeut = 0.0;
793     sigmaCharGlu  = 0.0;
794     sigmaNeutGlu  = comFacHat * 8.0 * alpEM * alpS / 9.0 / xW/(1-xW);
795   }
796
797 }
798
799 //--------------------------------------------------------------------------
800
801 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence. 
802
803 double Sigma2qq2squarksquark::sigmaHat() {
804
805   // In-pair must be same-sign
806   if (id1 * id2 < 0) return 0.0;
807   
808   // Check correct charge sum
809   if (isUD && abs(id1) %2 == abs(id2) % 2) return 0.0;
810   if (!isUD && abs(id1) % 2 != abs(id2) % 2) return 0.0;
811   if (!isUD && abs(id1) % 2 != abs(id3Sav) % 2) return 0.0;
812
813   // Coded sigma is for ud -> ~q~q'. Swap t and u for du -> ~q~q'.
814   swapTU = (isUD and abs(id1) % 2 == 0); 
815   int    idIn1A = (swapTU) ? abs(id2) : abs(id1);
816   int    idIn2A = (swapTU) ? abs(id1) : abs(id2);
817
818   // Auxiliary factors for use below
819   tGlu     = tH - m2Glu;
820   uGlu     = uH - m2Glu;
821   for (int i=1; i<= nNeut; i++) {
822     tNeut[i] = tH - m2Neut[i];
823     uNeut[i] = uH - m2Neut[i];
824     if (isUD && i <= 2) {
825       tChar[i] = tH - m2Char[i];
826       uChar[i] = uH - m2Char[i];
827     }
828   }
829
830   // Generation indices of incoming particles
831   int iGen1 = (abs(idIn1A)+1)/2;
832   int iGen2 = (abs(idIn2A)+1)/2;
833
834   // Initial values for pieces used for color-flow selection below
835   sumCt = 0.0;
836   sumCu = 0.0;
837   sumNt = 0.0;
838   sumNu = 0.0;
839   sumGt = 0.0;
840   sumGu = 0.0;
841   sumInterference = 0.0;
842
843   // Common factor for LR and RL contributions
844   double facTU =  uH*tH-s3*s4;
845
846   // Case A) Opposite-isospin: qq' -> ~d~u 
847   if ( isUD ) {
848
849     // t-channel Charginos
850     for (int k=1;k<=2;k++) {
851
852       // Skip if only including gluinos
853       if (settingsPtr->flag("SUSY:qq2squarksquark:onlyQCD")) continue;
854       
855       for (int l=1;l<=2;l++) {
856
857         // kl-dependent factor for LL and RR contributions
858         double facMS = sH*sqrt(m2Char[k]*m2Char[l]);
859
860         // Note: Ckl defined as in [Boz07] with sigmaChar factored out
861         // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
862         complex Ckl[3][3];
863         Ckl[1][1] = facMS
864           * coupSUSYPtr->LsudX[iGen4][iGen2][k]
865           * conj(coupSUSYPtr->LsduX[iGen3][iGen1][k])
866           * conj(coupSUSYPtr->LsudX[iGen4][iGen2][l])
867           * coupSUSYPtr->LsduX[iGen3][iGen1][l];
868         Ckl[1][2] = facTU
869           * coupSUSYPtr->RsudX[iGen4][iGen2][k]
870           * conj(coupSUSYPtr->LsduX[iGen3][iGen1][k])
871           * conj(coupSUSYPtr->RsudX[iGen4][iGen2][l])
872           * coupSUSYPtr->LsduX[iGen3][iGen1][l];
873         Ckl[2][1] = facTU 
874           * coupSUSYPtr->LsudX[iGen4][iGen2][k]
875           * conj(coupSUSYPtr->RsduX[iGen3][iGen1][k])
876           * conj(coupSUSYPtr->LsudX[iGen4][iGen2][l])
877           * coupSUSYPtr->RsduX[iGen3][iGen1][l];
878         Ckl[2][2] = facMS  
879           * coupSUSYPtr->RsudX[iGen4][iGen2][k]
880           * conj(coupSUSYPtr->RsduX[iGen3][iGen1][k])
881           * conj(coupSUSYPtr->RsudX[iGen4][iGen2][l])
882           * coupSUSYPtr->RsduX[iGen3][iGen1][l];
883
884         // Add to sum of t-channel charginos
885         sumCt += sigmaChar * real(Ckl[1][1] + Ckl[1][2] + Ckl[2][1] 
886                + Ckl[2][2]) / tChar[k] / tChar[l];
887
888       }
889     }
890     
891     // u-channel Neutralinos
892     for (int k=1;k<=nNeut;k++) {
893
894       // Skip if only including gluinos
895       if (settingsPtr->flag("SUSY:qq2squarksquark:onlyQCD")) continue;
896       
897       for (int l=1;l<=nNeut;l++) {
898
899         // kl-dependent factor for LL, RR contributions
900         double facMS = sH*sqrt(m2Neut[k]*m2Neut[l]);
901
902         // Note: Nkl defined as in [Boz07] with sigmaNeut factored out
903         // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
904         complex Nkl[3][3];
905         Nkl[1][1] = facMS
906           * conj(coupSUSYPtr->LsuuX[iGen4][iGen1][k])
907           * conj(coupSUSYPtr->LsddX[iGen3][iGen2][k])
908           * coupSUSYPtr->LsuuX[iGen4][iGen1][l]
909           * coupSUSYPtr->LsddX[iGen3][iGen2][l];
910         Nkl[1][2] = facTU 
911           * conj(coupSUSYPtr->LsuuX[iGen4][iGen1][k])
912           * conj(coupSUSYPtr->RsddX[iGen3][iGen2][k])
913           * coupSUSYPtr->LsuuX[iGen4][iGen1][l]
914           * coupSUSYPtr->RsddX[iGen3][iGen2][l];
915         Nkl[2][1] =  facTU 
916           * conj(coupSUSYPtr->RsuuX[iGen4][iGen1][k])
917           * conj(coupSUSYPtr->LsddX[iGen3][iGen2][k])
918           * coupSUSYPtr->RsuuX[iGen4][iGen1][l]
919           * coupSUSYPtr->LsddX[iGen3][iGen2][l];
920         Nkl[2][2] =  facMS 
921           * conj(coupSUSYPtr->RsuuX[iGen4][iGen1][k])
922           * conj(coupSUSYPtr->RsddX[iGen3][iGen2][k])
923           * coupSUSYPtr->RsuuX[iGen4][iGen1][l]
924           * coupSUSYPtr->RsddX[iGen3][iGen2][l];
925
926         // Add to sum of u-channel neutralinos  
927         sumNu += sigmaNeut / uNeut[k] / uNeut[l] 
928           * real(Nkl[1][1] + Nkl[1][2] + Nkl[2][1] + Nkl[2][2]);
929                 
930       }
931     }
932
933     // u-channel gluino
934     // Note: Gij defined as in [Boz07] with sigmaGlu factored out
935     // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
936     double Gij[3][3];
937     Gij[1][1] = norm(coupSUSYPtr->LsuuG[iGen4][iGen1]
938                 * coupSUSYPtr->LsddG[iGen3][iGen2]);
939     Gij[1][2] = norm(coupSUSYPtr->LsuuG[iGen4][iGen1]
940                 * coupSUSYPtr->RsddG[iGen3][iGen2]);
941     Gij[2][1] = norm(coupSUSYPtr->RsuuG[iGen4][iGen1]
942                 * coupSUSYPtr->LsddG[iGen3][iGen2]);
943     Gij[2][2] = norm(coupSUSYPtr->RsuuG[iGen4][iGen1]
944                 * coupSUSYPtr->RsddG[iGen3][iGen2]);
945     Gij[1][1] *= sH*m2Glu;
946     Gij[1][2] *= facTU;
947     Gij[2][1] *= facTU;
948     Gij[2][2] *= sH*m2Glu;
949     // Sum over polarizations
950     sumGu += sigmaGlu * (Gij[1][1] + Gij[1][2] + Gij[2][1] + Gij[2][2]) 
951       / pow2(uGlu); 
952
953     // chargino-neutralino interference
954     for (int k=1;k<=2;k++) {
955       for (int l=1;l<=nNeut;l++) {
956
957         // Skip if only including gluinos
958         if (settingsPtr->flag("SUSY:qq2squarksquark:onlyQCD")) continue;
959
960         // Note: CNkl defined as in [Boz07] with pi/sH2 factored out
961         // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
962         double CNkl[3][3];
963         CNkl[1][1] = real(coupSUSYPtr->LsudX[iGen4][iGen2][k] 
964                           * conj(coupSUSYPtr->LsduX[iGen3][iGen1][k]) 
965                           * coupSUSYPtr->LsuuX[iGen4][iGen1][l] 
966                           * coupSUSYPtr->LsddX[iGen3][iGen2][l]);
967         CNkl[1][2] = real(coupSUSYPtr->RsudX[iGen4][iGen2][k] 
968                           * conj(coupSUSYPtr->LsduX[iGen3][iGen1][k]) 
969                           * coupSUSYPtr->LsuuX[iGen4][iGen1][l] 
970                           * coupSUSYPtr->RsddX[iGen3][iGen2][l]);
971         CNkl[2][1] = real(coupSUSYPtr->LsudX[iGen4][iGen2][k] 
972                           * conj(coupSUSYPtr->RsduX[iGen3][iGen1][k]) 
973                           * coupSUSYPtr->RsuuX[iGen4][iGen1][l] 
974                           * coupSUSYPtr->LsddX[iGen3][iGen2][l]);
975         CNkl[2][2] = real(coupSUSYPtr->RsudX[iGen4][iGen2][k] 
976                           * conj(coupSUSYPtr->RsduX[iGen3][iGen1][k]) 
977                           * coupSUSYPtr->RsuuX[iGen4][iGen1][l] 
978                           * coupSUSYPtr->RsddX[iGen3][iGen2][l]);
979         CNkl[1][1] *= sH*sqrt(m2Char[k]*m2Neut[l]);
980         CNkl[1][2] *= uH*tH-s3*s4;
981         CNkl[2][1] *= uH*tH-s3*s4;
982         CNkl[2][2] *= sH*sqrt(m2Char[k]*m2Neut[l]);     
983         // Sum over polarizations
984         sumInterference += sigmaCharNeut * (CNkl[1][1] + CNkl[1][2] 
985                          + CNkl[2][1] + CNkl[2][2]) / tChar[k] / uNeut[l];
986       }
987     }
988
989     // chargino-gluino interference
990     for (int k=1;k<=2;k++) {
991
992       // Skip if only including gluinos
993       if (settingsPtr->flag("SUSY:qq2squarksquark:onlyQCD")) continue;
994         
995       // Note: CGk defined as in [Boz07] with sigmaCharGlu factored out
996       // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
997       double CGk[3][3];
998       CGk[1][1] = real(coupSUSYPtr->LsudX[iGen4][iGen2][k]
999                        * conj(coupSUSYPtr->LsduX[iGen3][iGen1][k])
1000                        * conj(coupSUSYPtr->LsuuG[iGen4][iGen1])
1001                        * conj(coupSUSYPtr->LsddG[iGen3][iGen2]));
1002       CGk[1][2] = real(coupSUSYPtr->RsudX[iGen4][iGen2][k]
1003                        * conj(coupSUSYPtr->LsduX[iGen3][iGen1][k])
1004                        * conj(coupSUSYPtr->LsuuG[iGen4][iGen1])
1005                        * conj(coupSUSYPtr->RsddG[iGen3][iGen2]));
1006       CGk[2][1] = real(coupSUSYPtr->LsudX[iGen4][iGen2][k]
1007                        * conj(coupSUSYPtr->RsduX[iGen3][iGen1][k])
1008                        * conj(coupSUSYPtr->RsuuG[iGen4][iGen1])
1009                        * conj(coupSUSYPtr->LsddG[iGen3][iGen2]));
1010       CGk[2][2] = real(coupSUSYPtr->RsudX[iGen4][iGen2][k]
1011                        * conj(coupSUSYPtr->RsduX[iGen3][iGen1][k])
1012                        * conj(coupSUSYPtr->RsuuG[iGen4][iGen1])
1013                        * conj(coupSUSYPtr->RsddG[iGen3][iGen2]));
1014       CGk[1][1] *= sH*sqrt(m2Glu*m2Char[k]);
1015       CGk[1][2] *= uH*tH-s3*s4;
1016       CGk[2][1] *= uH*tH-s3*s4;
1017       CGk[2][2] *= sH*sqrt(m2Glu*m2Char[k]);
1018       // Sum over polarizations
1019       sumInterference += sigmaGlu * (CGk[1][1] + CGk[1][2] + CGk[2][1] 
1020         + CGk[2][2]) / uGlu / tChar[k]; 
1021     }    
1022   }
1023
1024   // Case B) Same-isospin: qq' -> ~d~d , ~u~u
1025   else {    
1026
1027     // t-channel + u-channel Neutralinos + t/u interference
1028     for (int k=1;k<=nNeut;k++) {
1029
1030       // Skip if only including gluinos
1031       if (settingsPtr->flag("SUSY:qq2squarksquark:onlyQCD")) continue;
1032
1033       for (int l=1;l<=nNeut;l++) {
1034
1035         // kl-dependent factor for LL and RR contributions
1036         double facMS = sH * particleDataPtr->m0(coupSUSYPtr->idNeut(k))
1037           * particleDataPtr->m0(coupSUSYPtr->idNeut(l));
1038
1039         // Note: Nxkl defined as in [Boz07] with sigmaNeut factored out
1040         // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
1041         complex NTkl[3][3], NUkl[3][3], NTUkl[3][3];
1042         NTkl[1][1] = facMS 
1043           * conj(coupSUSYPtr->getLsqqX(iGen4,idIn2A,k))
1044           * conj(coupSUSYPtr->getLsqqX(iGen3,idIn1A,k))
1045           * coupSUSYPtr->getLsqqX(iGen4,idIn2A,l)
1046           * coupSUSYPtr->getLsqqX(iGen3,idIn1A,l);
1047         NTkl[1][2] = facTU 
1048           * conj(coupSUSYPtr->getRsqqX(iGen4,idIn2A,k))
1049           * conj(coupSUSYPtr->getLsqqX(iGen3,idIn1A,k))
1050           * coupSUSYPtr->getRsqqX(iGen4,idIn2A,l)
1051           * coupSUSYPtr->getLsqqX(iGen3,idIn1A,l);
1052         NTkl[2][1] = facTU 
1053           * conj(coupSUSYPtr->getLsqqX(iGen4,idIn2A,k))
1054           * conj(coupSUSYPtr->getRsqqX(iGen3,idIn1A,k))
1055           * coupSUSYPtr->getLsqqX(iGen4,idIn2A,l)
1056           * coupSUSYPtr->getRsqqX(iGen3,idIn1A,l);
1057         NTkl[2][2] = facMS  
1058           * conj(coupSUSYPtr->getRsqqX(iGen4,idIn2A,k))
1059           * conj(coupSUSYPtr->getRsqqX(iGen3,idIn1A,k))
1060           * coupSUSYPtr->getRsqqX(iGen4,idIn2A,l)
1061           * coupSUSYPtr->getRsqqX(iGen3,idIn1A,l);
1062         NUkl[1][1] = facMS
1063           * conj(coupSUSYPtr->getLsqqX(iGen3,idIn2A,k))
1064           * conj(coupSUSYPtr->getLsqqX(iGen4,idIn1A,k))
1065           * coupSUSYPtr->getLsqqX(iGen3,idIn2A,l)
1066           * coupSUSYPtr->getLsqqX(iGen4,idIn1A,l);
1067         NUkl[1][2] = facTU 
1068           * conj(coupSUSYPtr->getRsqqX(iGen3,idIn2A,k))
1069           * conj(coupSUSYPtr->getLsqqX(iGen4,idIn1A,k))
1070           * coupSUSYPtr->getRsqqX(iGen3,idIn2A,l)
1071           * coupSUSYPtr->getLsqqX(iGen4,idIn1A,l);
1072         NUkl[2][1] = facTU
1073           * conj(coupSUSYPtr->getLsqqX(iGen3,idIn2A,k))
1074           * conj(coupSUSYPtr->getRsqqX(iGen4,idIn1A,k))
1075           * coupSUSYPtr->getLsqqX(iGen3,idIn2A,l)
1076           * coupSUSYPtr->getRsqqX(iGen4,idIn1A,l);
1077         NUkl[2][2] = facMS
1078           * conj(coupSUSYPtr->getRsqqX(iGen3,idIn2A,k))
1079           * conj(coupSUSYPtr->getRsqqX(iGen4,idIn1A,k))
1080           * coupSUSYPtr->getRsqqX(iGen3,idIn2A,l)
1081           * coupSUSYPtr->getRsqqX(iGen4,idIn1A,l);
1082         NTUkl[1][1] = facMS  
1083           * real( conj(coupSUSYPtr->getLsqqX(iGen4,idIn2A,k))
1084                   * conj(coupSUSYPtr->getLsqqX(iGen3,idIn1A,k))
1085                   * coupSUSYPtr->getLsqqX(iGen3,idIn2A,l)
1086                   * coupSUSYPtr->getLsqqX(iGen4,idIn1A,l) );
1087         NTUkl[1][2] = facTU 
1088           * real( conj(coupSUSYPtr->getRsqqX(iGen4,idIn2A,k))
1089                   * conj(coupSUSYPtr->getLsqqX(iGen3,idIn1A,k))
1090                   * coupSUSYPtr->getRsqqX(iGen3,idIn2A,l)
1091                   * coupSUSYPtr->getLsqqX(iGen4,idIn1A,l) );
1092         NTUkl[2][1] = facTU 
1093           * real( conj(coupSUSYPtr->getLsqqX(iGen4,idIn2A,k))
1094                   * conj(coupSUSYPtr->getRsqqX(iGen3,idIn1A,k))
1095                   * coupSUSYPtr->getLsqqX(iGen3,idIn2A,l)
1096                   * coupSUSYPtr->getRsqqX(iGen4,idIn1A,l) );
1097         NTUkl[2][2] = facMS  
1098           * real( conj(coupSUSYPtr->getRsqqX(iGen4,idIn2A,k))
1099                   * conj(coupSUSYPtr->getRsqqX(iGen3,idIn1A,k))
1100                   * coupSUSYPtr->getRsqqX(iGen3,idIn2A,l)
1101                   * coupSUSYPtr->getRsqqX(iGen4,idIn1A,l) );
1102           
1103         // Add to sums
1104         sumNt += sigmaNeut / tNeut[k] / tNeut[l] 
1105           * real(NTkl[1][1] + NTkl[1][2] + NTkl[2][1] + NTkl[2][2]);
1106         sumNu += sigmaNeut / uNeut[k] / uNeut[l]
1107           * real(NUkl[1][1] + NUkl[1][2] + NUkl[2][1] + NUkl[2][2]);
1108         sumInterference += 2.0 / 3.0 * sigmaNeut 
1109           * real(NTUkl[1][1] + NTUkl[1][2] + NTUkl[2][1] + NTUkl[2][2])
1110           / tNeut[k] / uNeut[l];
1111       }     
1112
1113       // Neutralino / Gluino interference
1114
1115       // k-dependent factor for LL and RR contributions
1116       double facMS = sH * particleDataPtr->m0(coupSUSYPtr->idNeut(k))
1117         * particleDataPtr->m0(1000021);
1118       
1119       // Note: Nxkl defined as in [Boz07] with sigmaNeutGlu factored out
1120       // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
1121       complex NGA[3][3], NGB[3][3];
1122       NGA[1][1] = facMS
1123         * real( conj(coupSUSYPtr->getLsqqX(iGen4,idIn2A,k))
1124                 * conj(coupSUSYPtr->getLsqqX(iGen3,idIn1A,k))
1125                 * conj(coupSUSYPtr->getLsqqG(iGen3,idIn2A))
1126                 * conj(coupSUSYPtr->getLsqqG(iGen4,idIn1A)) );
1127       NGA[1][2] = facTU
1128         * real( conj(coupSUSYPtr->getRsqqX(iGen4,idIn2A,k))
1129                 * conj(coupSUSYPtr->getLsqqX(iGen3,idIn1A,k))
1130                 * conj(coupSUSYPtr->getLsqqG(iGen3,idIn2A))
1131                 * conj(coupSUSYPtr->getRsqqG(iGen4,idIn1A)) );
1132       NGA[2][1] = facTU
1133         * real( conj(coupSUSYPtr->getLsqqX(iGen4,idIn2A,k))
1134                 * conj(coupSUSYPtr->getRsqqX(iGen3,idIn1A,k))
1135                 * conj(coupSUSYPtr->getRsqqG(iGen3,idIn2A))
1136                 * conj(coupSUSYPtr->getLsqqG(iGen4,idIn1A)) );
1137       NGA[2][2] = facMS
1138         * real( conj(coupSUSYPtr->getRsqqX(iGen4,idIn2A,k))
1139                 * conj(coupSUSYPtr->getRsqqX(iGen3,idIn1A,k))
1140                 * conj(coupSUSYPtr->getRsqqG(iGen3,idIn2A))
1141                 * conj(coupSUSYPtr->getRsqqG(iGen4,idIn1A)) );
1142       NGB[1][1] = facMS
1143         * real( conj(coupSUSYPtr->getLsqqX(iGen3,idIn2A,k))
1144                 * conj(coupSUSYPtr->getLsqqX(iGen4,idIn1A,k))
1145                 * conj(coupSUSYPtr->getLsqqG(iGen4,idIn2A))
1146                 * conj(coupSUSYPtr->getLsqqG(iGen3,idIn1A)) );
1147       NGB[1][2] = facMS
1148         * real( conj(coupSUSYPtr->getRsqqX(iGen3,idIn2A,k))
1149                 * conj(coupSUSYPtr->getLsqqX(iGen4,idIn1A,k))
1150                 * conj(coupSUSYPtr->getRsqqG(iGen4,idIn2A))
1151                 * conj(coupSUSYPtr->getLsqqG(iGen3,idIn1A)) );
1152       NGB[2][1] = facMS
1153         * real( conj(coupSUSYPtr->getLsqqX(iGen3,idIn2A,k))
1154                 * conj(coupSUSYPtr->getRsqqX(iGen4,idIn1A,k))
1155                 * conj(coupSUSYPtr->getLsqqG(iGen4,idIn2A))
1156                 * conj(coupSUSYPtr->getRsqqG(iGen3,idIn1A)) );
1157       NGB[2][2] = facMS
1158         * real( conj(coupSUSYPtr->getRsqqX(iGen3,idIn2A,k))
1159                 * conj(coupSUSYPtr->getRsqqX(iGen4,idIn1A,k))
1160                 * conj(coupSUSYPtr->getRsqqG(iGen4,idIn2A))
1161                 * conj(coupSUSYPtr->getRsqqG(iGen3,idIn1A)) );
1162
1163       // Add to sums
1164       sumInterference += sigmaNeutGlu * 
1165         ( real(NGA[1][1] + NGA[1][2] + NGA[2][1] + NGA[2][2]) 
1166         / tNeut[k] / uGlu
1167         + real(NGB[1][1] + NGB[1][2] + NGB[2][1] + NGB[2][2]) 
1168         / uNeut[k] / tGlu );
1169     }
1170     
1171     // t-channel + u-channel Gluinos + t/u interference
1172
1173     // factor for LL and RR contributions
1174     double facMS = sH * m2Glu;
1175
1176     // Note: GT, GU defined as in [Boz07] with sigmaGlu factored out
1177     // [1][1] = LL, [1][2] = LR, [2][1] = RL, [2][2] = RR
1178     complex GT[3][3], GU[3][3], GTU[3][3];
1179     GT[1][1] = facMS 
1180       * norm(coupSUSYPtr->getLsqqG(iGen4,idIn2A) 
1181       * coupSUSYPtr->getLsqqG(iGen3,idIn1A));
1182     GT[1][2] = facTU
1183       * norm(coupSUSYPtr->getRsqqG(iGen4,idIn2A) 
1184       * coupSUSYPtr->getLsqqG(iGen3,idIn1A));
1185     GT[2][1] = facTU
1186       * norm(coupSUSYPtr->getLsqqG(iGen4,idIn2A) 
1187       * coupSUSYPtr->getRsqqG(iGen3,idIn1A));
1188     GT[2][2] = facMS
1189       * norm(coupSUSYPtr->getRsqqG(iGen4,idIn2A) 
1190       * coupSUSYPtr->getRsqqG(iGen3,idIn1A));
1191     GU[1][1] = facMS 
1192       * norm(coupSUSYPtr->getLsqqG(iGen3,idIn2A) 
1193       * coupSUSYPtr->getLsqqG(iGen4,idIn1A));
1194     GU[1][2] = facTU
1195       * norm(coupSUSYPtr->getLsqqG(iGen3,idIn2A) 
1196       * coupSUSYPtr->getRsqqG(iGen4,idIn1A));
1197     GU[2][1] = facTU
1198       * norm(coupSUSYPtr->getRsqqG(iGen3,idIn2A) 
1199       * coupSUSYPtr->getLsqqG(iGen4,idIn1A));
1200     GU[2][2] = facMS
1201       * norm(coupSUSYPtr->getRsqqG(iGen3,idIn2A) 
1202       * coupSUSYPtr->getRsqqG(iGen4,idIn1A));
1203
1204     GTU[1][1] = facMS 
1205       * real(coupSUSYPtr->getLsqqG(iGen3,idIn1A) 
1206              * coupSUSYPtr->getLsqqG(iGen4,idIn2A)
1207              * conj(coupSUSYPtr->getLsqqG(iGen3,idIn2A))
1208              * conj(coupSUSYPtr->getLsqqG(iGen4,idIn1A)) );
1209
1210     GTU[1][2] = facTU 
1211       * real(coupSUSYPtr->getLsqqG(iGen3,idIn1A) 
1212              * coupSUSYPtr->getRsqqG(iGen4,idIn2A)
1213              * conj(coupSUSYPtr->getRsqqG(iGen3,idIn2A))
1214              * conj(coupSUSYPtr->getLsqqG(iGen4,idIn1A)) );
1215
1216     GTU[2][1] = facTU 
1217       * real(coupSUSYPtr->getRsqqG(iGen3,idIn1A) 
1218              * coupSUSYPtr->getLsqqG(iGen4,idIn2A)
1219              * conj(coupSUSYPtr->getLsqqG(iGen3,idIn2A))
1220              * conj(coupSUSYPtr->getRsqqG(iGen4,idIn1A)) );
1221
1222     GTU[2][2] = facMS 
1223       * real(coupSUSYPtr->getRsqqG(iGen3,idIn1A) 
1224              * coupSUSYPtr->getRsqqG(iGen4,idIn2A)
1225              * conj(coupSUSYPtr->getRsqqG(iGen3,idIn2A))
1226              * conj(coupSUSYPtr->getRsqqG(iGen4,idIn1A)) );
1227  
1228     // Add to sums
1229     sumGt += sigmaGlu * real(GT[1][1] + GT[1][2] + GT[2][1] + GT[2][2])
1230       / pow2(tGlu) ;
1231     sumGu += sigmaGlu * real(GU[1][1] + GU[1][2] + GU[2][1] + GU[2][2])
1232       / pow2(uGlu) ;
1233     sumInterference += - 2.0 / 3.0 * sigmaGlu
1234       * real(GTU[1][1] + GTU[1][2] + GTU[2][1] + GTU[2][2])
1235       / tGlu / uGlu;    
1236
1237   }
1238
1239   // Cross section
1240   double sigma = sumNt + sumNu + sumCt + sumCu + sumGt + sumGu 
1241     + sumInterference;
1242
1243   // Identical particles?
1244   if (id3Sav == id4Sav) sigma /= 2.0;
1245
1246   // Return answer.
1247   return sigma;
1248
1249 }
1250
1251 //--------------------------------------------------------------------------
1252
1253 // Select identity, colour and anticolour.
1254
1255 void Sigma2qq2squarksquark::setIdColAcol() {
1256
1257   // Set flavours.
1258   if (id1 > 0 && id2 > 0) {
1259     setId( id1, id2, id3Sav, id4Sav);
1260   } else {
1261     // 1,2 -> -3,-4
1262     setId( id1, id2,-id3Sav,-id4Sav);
1263   }
1264
1265   // Coded sigma is for ud -> ~q~q'. Swap t and u for du -> ~q~q'.
1266   swapTU = (isUD and abs(id1) % 2 == 0); 
1267
1268   // Select colour flow topology 
1269   // A: t-channel neutralino, t-channel chargino, or u-channel gluino
1270   double fracA = sumNt + sumCt + sumGu 
1271     / (sumNt + sumNu + sumCt + sumCu + sumGt + sumGu);
1272   if (swapTU) fracA = 1.0 - fracA;
1273   setColAcol( 1, 0, 2, 0, 1, 0, 2, 0);
1274   // B: t-channel gluino or u-channel neutralino 
1275   if (rndmPtr->flat() > fracA) setColAcol( 1, 0, 2, 0, 2, 0, 1, 0);
1276
1277   // Switch to anti-colors if antiquarks
1278   if (id1 < 0 || id2 < 0) swapColAcol();
1279
1280 }
1281
1282 //==========================================================================
1283
1284 // Sigma2qqbar2squarkantisquark
1285 // Cross section for qqbar-initiated squark-antisquark production
1286
1287 //--------------------------------------------------------------------------
1288
1289 // Initialize process. 
1290   
1291 void Sigma2qqbar2squarkantisquark::initProc() {
1292
1293   //Typecast to the correct couplings
1294   coupSUSYPtr = (CoupSUSY*) couplingsPtr;
1295
1296   // Extract isospin and mass-ordering indices
1297   iGen3 = 3*(abs(id3Sav)/2000000) + (abs(id3Sav)%10+1)/2;
1298   iGen4 = 3*(abs(id4Sav)/2000000) + (abs(id4Sav)%10+1)/2;
1299
1300   // Is this a ~u_i ~d*_j, ~d_i ~u*_j final state or ~d_i ~d*_j, ~u_i ~u*_j
1301   if (abs(id3Sav) % 2 == abs(id4Sav) % 2) isUD = false;
1302   else isUD = true;
1303
1304   // Derive name
1305   nameSave = "q qbar' -> "+particleDataPtr->name(abs(id3Sav))+" "+
1306     particleDataPtr->name(-abs(id4Sav));
1307   if (isUD && abs(id3Sav) != abs(id4Sav)) nameSave +=" + c.c.";
1308
1309   // Count 5 neutralinos in NMSSM
1310   nNeut = (coupSUSYPtr->isNMSSM ? 5 : 4);
1311
1312   // Store mass squares of all possible internal propagator lines
1313   m2Glu     = pow2(particleDataPtr->m0(1000021));
1314   m2Neut.resize(nNeut+1);
1315   for (int iNeut=1;iNeut<=nNeut;iNeut++) 
1316     m2Neut[iNeut] = pow2(particleDataPtr->m0(coupSUSYPtr->idNeut(iNeut)));
1317   
1318   // Set sizes of some arrays to be used below
1319   tNeut.resize(nNeut+1);
1320   uNeut.resize(nNeut+1);
1321
1322   // Shorthand for Weak mixing
1323   xW = coupSUSYPtr->sin2W;
1324
1325   // Secondary open width fraction.
1326   openFracPair = particleDataPtr->resOpenFrac(id3Sav, id4Sav);
1327
1328 }
1329
1330 //--------------------------------------------------------------------------
1331
1332 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour. 
1333
1334 void Sigma2qqbar2squarkantisquark::sigmaKin() {
1335
1336   // Z/W propagator
1337   if (! isUD) {
1338     double sV= sH - pow2(coupSUSYPtr->mZpole);
1339     double d = pow2(sV) + pow2(coupSUSYPtr->mZpole * coupSUSYPtr->wZpole);
1340     propZW   = complex( sV / d, coupSUSYPtr->mZpole * coupSUSYPtr->wZpole / d);
1341   } else {
1342     double sV= sH - pow2(coupSUSYPtr->mWpole);
1343     double d = pow2(sV) + pow2(coupSUSYPtr->mWpole * coupSUSYPtr->wWpole);
1344     propZW   = complex( sV / d, coupSUSYPtr->mWpole * coupSUSYPtr->wWpole / d);
1345   }
1346
1347   // Flavor-independent pre-factors
1348   double comFacHat = M_PI/sH2 * openFracPair;
1349
1350   sigmaEW       = comFacHat * pow2(alpEM);
1351   sigmaGlu      = comFacHat * 2.0 * pow2(alpS) / 9.0;
1352   sigmaEWG      = comFacHat * 8.0 * alpEM * alpS / 9.0;
1353
1354 }
1355
1356 //--------------------------------------------------------------------------
1357
1358 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence. 
1359
1360 double Sigma2qqbar2squarkantisquark::sigmaHat() {
1361
1362   // In-pair must be opposite-sign
1363   if (id1 * id2 > 0) return 0.0;
1364   
1365   // Check correct charge sum
1366   if (isUD && abs(id1) %2 == abs(id2) % 2) return 0.0;
1367   if (!isUD && abs(id1) % 2 != abs(id2) % 2) return 0.0;
1368
1369   // Check if using QCD diagrams only
1370   bool onlyQCD = settingsPtr->flag("SUSY:qqbar2squarkantisquark:onlyQCD");
1371
1372   // Coded UD sigma is for udbar -> ~u~d'*. Swap t<->u for dbaru -> ~u~d'*.
1373   swapTU = (isUD and abs(id1) % 2 != 0); 
1374
1375   // Coded QQ sigma is for qqbar -> ~q~q*. Swap t<->u for qbarq -> ~q~q*.
1376   if (!isUD && id1 < 0) swapTU = true;
1377
1378   // Generation indices of incoming particles
1379   int idIn1A = (swapTU) ? abs(id2) : abs(id1);
1380   int idIn2A = (swapTU) ? abs(id1) : abs(id2);
1381   int iGen1  = (idIn1A+1)/2;
1382   int iGen2  = (idIn2A+1)/2;  
1383
1384   // Auxiliary factors for use below  
1385   tGlu     = tH - m2Glu;
1386   uGlu     = uH - m2Glu;
1387   for (int i=1; i<= nNeut; i++) {
1388     tNeut[i] = tH - m2Neut[i];
1389     uNeut[i] = uH - m2Neut[i];
1390   }
1391
1392   // Initial values for pieces used for color-flow selection below
1393   sumColS   = 0.0;
1394   sumColT   = 0.0;
1395   sumInterference = 0.0;
1396
1397   // Common factor for LR and RL contributions
1398   double facTU =  uH*tH-s3*s4;
1399
1400   // Case A) Opposite-isospin: udbar -> ~u~d* 
1401   if ( isUD ) {
1402
1403     // s-channel W contribution (only contributes to LL helicities)
1404     if ( !onlyQCD ) {
1405       sumColS += sigmaEW / 16.0 / pow2(xW) / pow2(1.0-xW)
1406         * norm(conj(coupSUSYPtr->LudW[iGen1][iGen2])
1407                * coupSUSYPtr->LsusdW[iGen3][iGen4]) * facTU  
1408         * norm(propZW);
1409     }
1410
1411     // t-channel gluino contributions
1412     double GT[3][3];
1413     double facLR = m2Glu * sH;
1414     // LL, LR, RL, RR
1415     GT[1][1] = facTU * norm(conj(coupSUSYPtr->LsddG[iGen4][iGen2])
1416                             *coupSUSYPtr->LsuuG[iGen3][iGen1]);
1417     GT[1][2] = facLR * norm(conj(coupSUSYPtr->RsddG[iGen4][iGen2])
1418                             *coupSUSYPtr->LsuuG[iGen3][iGen1]);
1419     GT[2][1] = facLR * norm(conj(coupSUSYPtr->LsddG[iGen4][iGen2])
1420                             *coupSUSYPtr->RsuuG[iGen3][iGen1]);
1421     GT[2][2] = facTU * norm(conj(coupSUSYPtr->RsddG[iGen4][iGen2])
1422                             *coupSUSYPtr->RsuuG[iGen3][iGen1]);
1423     // leading color flow for t-channel gluino is annihilation-like
1424     sumColS += sigmaGlu / pow2(tGlu)
1425       * (GT[1][1] + GT[1][2] + GT[2][1] + GT[2][2]);
1426       
1427     // W-Gluino interference (only contributes to LL helicities)
1428     if ( !onlyQCD ) {
1429       sumColS += sigmaEWG / 4.0 / xW / (1-xW) 
1430         * real(conj(coupSUSYPtr->LsuuG[iGen3][iGen1])
1431                * coupSUSYPtr->LsddG[iGen4][iGen2]
1432                * conj(coupSUSYPtr->LudW[iGen1][iGen2])
1433                * coupSUSYPtr->LsusdW[iGen3][iGen4]) * facTU 
1434         / tGlu * sqrt(norm(propZW));
1435     }
1436
1437     // t-channel neutralinos
1438     // NOT YET IMPLEMENTED !
1439
1440   }
1441
1442   // Case B) Same-isospin: qqbar -> ~d~d* , ~u~u*
1443   else {    
1444     
1445     double eQ  = (idIn1A % 2 == 0) ? 2./3. : 1./3. ;
1446     double eSq = (abs(id3Sav) % 2 == 0) ? 2./3. : 1./3. ;
1447
1448     // s-channel gluon (strictly flavor-diagonal)
1449     if (abs(id3Sav) == abs(id4Sav) && abs(id1) == abs(id2)) {
1450       // Factor 2 since contributes to both ha != hb helicities
1451       sumColT += 2. * sigmaGlu * facTU / pow2(sH);      
1452     }
1453
1454     // t-channel gluino (only for in-isospin = out-isospin). 
1455     if (eQ == eSq) {
1456       // Sum over helicities.     
1457       double GT[3][3];
1458       double facLR = sH * m2Glu;
1459       GT[1][1] = facTU * norm(coupSUSYPtr->getLsqqG(iGen3,idIn1A)
1460                               * conj(coupSUSYPtr->getLsqqG(iGen4,idIn2A)));
1461       GT[1][2] = facLR * norm(coupSUSYPtr->getLsqqG(iGen3,idIn1A)
1462                               * conj(coupSUSYPtr->getRsqqG(iGen4,idIn2A)));
1463       GT[2][1] = facLR * norm(coupSUSYPtr->getRsqqG(iGen3,idIn1A)
1464                               * conj(coupSUSYPtr->getLsqqG(iGen4,idIn2A)));
1465       GT[2][2] = facTU * norm(coupSUSYPtr->getRsqqG(iGen3,idIn1A)
1466                               * conj(coupSUSYPtr->getRsqqG(iGen4,idIn2A)));    
1467       // Add contribution to color topology: S
1468       sumColS += sigmaGlu / pow2(tGlu)
1469         * ( GT[1][1] + GT[1][2] + GT[2][1] + GT[2][2]);
1470       
1471       // gluon-gluino interference (strictly flavor-diagonal)
1472       if (abs(id3Sav) == abs(id4Sav) && abs(id1) == abs(id2)) {
1473         double GG11, GG22;
1474         GG11 = - facTU * 2./3. 
1475              * real( conj(coupSUSYPtr->getLsqqG(iGen3,idIn1A))
1476              * coupSUSYPtr->getLsqqG(iGen4,idIn2A));
1477         GG22 = - facTU * 2./3. 
1478              * real( conj(coupSUSYPtr->getRsqqG(iGen3,idIn1A))
1479              * coupSUSYPtr->getRsqqG(iGen4,idIn2A)); 
1480         // Sum over two contributing helicities
1481         sumInterference += sigmaGlu / sH / tGlu
1482           * ( GG11 + GG22 );
1483       }
1484
1485     }
1486
1487     // Skip the rest if only including QCD diagrams
1488     if (onlyQCD) return sumColT+sumColS+sumInterference;
1489
1490     // s-channel photon (strictly flavor-diagonal) and Z/gamma interference
1491     if (abs(id3Sav) == abs(id4Sav) && abs(id1) == abs(id2)) {
1492
1493       // gamma
1494       // Factor 2 since contributes to both ha != hb helicities
1495       sumColS += 2. * pow2(eQ) * pow2(eSq) * sigmaEW * facTU / pow2(sH);
1496
1497       // Z/gamma interference
1498       double CsqZ = real(coupSUSYPtr->LsusuZ[iGen3][iGen4] 
1499                          + coupSUSYPtr->RsusuZ[iGen3][iGen4]);
1500       if (abs(id3Sav)%2 != 0) CsqZ = real(coupSUSYPtr->LsdsdZ[iGen3][iGen4] 
1501                                           + coupSUSYPtr->RsdsdZ[iGen3][iGen4]);
1502       sumColS += eQ * eSq * sigmaEW * facTU / 2.0 / xW / (1.-xW) 
1503         * sqrt(norm(propZW)) / sH * CsqZ
1504         * (coupSUSYPtr->LqqZ[idIn1A] + coupSUSYPtr->LqqZ[idIn2A]);
1505       
1506       // Gluino/gamma interference (only for same-isospin)
1507       if (eQ == eSq) {
1508         double CsqG11 = real(conj(coupSUSYPtr->LsuuG[iGen3][iGen1]) 
1509                              *coupSUSYPtr->LsuuG[iGen4][iGen2]);
1510         double CsqG22 = real(conj(coupSUSYPtr->RsuuG[iGen3][iGen1]) 
1511                              *coupSUSYPtr->RsuuG[iGen4][iGen2]);
1512         if (id3Sav%2 != 0) {
1513           CsqG11 = real(conj(coupSUSYPtr->LsddG[iGen3][iGen1]) 
1514                         *coupSUSYPtr->LsddG[iGen4][iGen2]);
1515           CsqG22 = real(conj(coupSUSYPtr->RsddG[iGen3][iGen1]) 
1516                         *coupSUSYPtr->RsddG[iGen4][iGen2]);
1517         }
1518         sumColS += eQ * eSq * sigmaEWG * facTU
1519           * (CsqG11 + CsqG22) / sH / tGlu; 
1520       }
1521     }
1522     
1523     // s-channel Z (only for q flavor = qbar flavor)
1524     if (abs(id1) == abs(id2)) {
1525       double CsqZ = norm(coupSUSYPtr->LsusuZ[iGen3][iGen4] 
1526                          + coupSUSYPtr->RsusuZ[iGen3][iGen4]);
1527       if (abs(id3Sav)%2 != 0) CsqZ = norm(coupSUSYPtr->LsdsdZ[iGen3][iGen4] 
1528                                           + coupSUSYPtr->RsdsdZ[iGen3][iGen4]);
1529       sumColS += sigmaEW * facTU / 16.0 / pow2(xW) / pow2(1.0-xW)
1530         * norm(propZW) * CsqZ * ( pow2(coupSUSYPtr->LqqZ[idIn1A]) 
1531         + pow2(coupSUSYPtr->RqqZ[idIn1A]) );
1532
1533       // Z/gluino interference (only for in-isospin = out-isospin)
1534       if (eQ == eSq) {
1535         double GZ11 = real(conj(coupSUSYPtr->getLsqqG(iGen3,idIn1A))
1536                            *coupSUSYPtr->getLsqqG(iGen4,idIn2A)          
1537                            *(coupSUSYPtr->getLsqsqZ(id3Sav,id4Sav)
1538                              +coupSUSYPtr->getRsqsqZ(id3Sav,id4Sav)))
1539           *coupSUSYPtr->LqqZ[idIn1A];
1540         double GZ22 = real(conj(coupSUSYPtr->getRsqqG(iGen3,idIn1A))
1541                            *coupSUSYPtr->getRsqqG(iGen4,idIn2A)
1542                            *(coupSUSYPtr->getLsqsqZ(id3Sav,id4Sav)
1543                              +coupSUSYPtr->getRsqsqZ(id3Sav,id4Sav)))
1544           *coupSUSYPtr->RqqZ[idIn1A];
1545         sumColS += sigmaEWG * facTU / 4.0 / xW / (1.-xW) 
1546           * ( GZ11 + GZ22 ) * sqrt(norm(propZW)) / tGlu;        
1547       }
1548     }
1549     
1550     // t-channel neutralinos
1551     // NOT YET IMPLEMENTED !
1552     
1553   }
1554
1555   // Cross section
1556   double sigma = sumColS + sumColT + sumInterference;
1557
1558   // Return answer.
1559   return sigma;
1560   
1561 }
1562
1563 //--------------------------------------------------------------------------
1564
1565 // Select identity, colour and anticolour.
1566
1567 void Sigma2qqbar2squarkantisquark::setIdColAcol() {
1568
1569   // Check if charge conjugate final state?
1570   isCC = false;
1571   if (isUD && ( (id1-1)%2 < 0 || (id2-1)%2 < 0 )) isCC = true;
1572   
1573   //check if charge conjugate
1574   id3 = (isCC) ? -id3Sav : id3Sav;
1575   id4 = (isCC) ? -id4Sav : id4Sav;
1576
1577   // Set flavours.
1578   setId( id1, id2, id3, id4);                 
1579
1580   // Coded UD sigma is for udbar -> ~u~d'*. Swap t<->u for dbaru -> ~u~d'*.
1581   // Coded QQ sigma is for qqbar -> ~q~q*. Swap t<->u for qbarq -> ~q~q*.
1582   if (isUD) {
1583     swapTU = (abs(id1) % 2 != 0); 
1584   } else {
1585     swapTU = (id1 < 0);
1586   }
1587
1588   // Select colour flow topology 
1589   double R = rndmPtr->flat();
1590   double fracS = sumColS / (sumColS + sumColT) ;
1591   // S: color flow as in S-channel singlet
1592   if (R < fracS) {
1593     setColAcol( 1, 0, 0, 1, 2, 0, 0, 2);
1594     if (swapTU) setColAcol( 0, 1, 1, 0, 2, 0, 0, 2);
1595   } 
1596   // T: color flow as in T-channel singlet
1597   else {
1598     setColAcol( 1, 0, 0, 2, 1, 0, 0, 2);
1599     if (swapTU) setColAcol( 0, 1, 2, 0, 2, 0, 0, 1);
1600   }
1601
1602   if (isCC) swapColAcol();
1603
1604 }
1605
1606 //==========================================================================
1607
1608 // Sigma2gg2squarkantisquark
1609 // Cross section for gg-initiated squark-antisquark production
1610
1611 //--------------------------------------------------------------------------
1612
1613 // Initialize process. 
1614   
1615 void Sigma2gg2squarkantisquark::initProc() {
1616
1617   //Typecast to the correct couplings
1618   coupSUSYPtr = (CoupSUSY*) couplingsPtr;
1619
1620   // Process Name
1621   nameSave = "g g -> "+particleDataPtr->name(abs(id3Sav))+" "
1622     +particleDataPtr->name(-abs(id4Sav)); 
1623
1624   // Squark pole mass
1625   m2Sq = pow2(particleDataPtr->m0(id3Sav));
1626   
1627   // Secondary open width fraction.
1628   openFracPair = particleDataPtr->resOpenFrac(id3Sav, id4Sav);
1629
1630 }
1631
1632 //--------------------------------------------------------------------------
1633
1634 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour. 
1635
1636 void Sigma2gg2squarkantisquark::sigmaKin() {
1637
1638   // Modified Mandelstam variables for massive kinematics with m3 = m4.
1639   // tHSq = tHat - m_squark^2; uHSq = uHat - m_squark^2. 
1640   //  double s34Avg = 0.5 * (s3 + s4) - 0.25 * pow2(s3 - s4) / sH; 
1641   double tHSq    = -0.5 * (sH - tH + uH);
1642   double uHSq    = -0.5 * (sH + tH - uH); 
1643   // ! (NEED TO CHECK THAT THESE APPLIED CORRECTLY BELOW)   ! 
1644   // ! (PRELIMINARY CROSS-CHECKS WITH PYTHIA 6 COME OUT OK) !
1645
1646   // Helicity-independent prefactor
1647   double comFacHat = M_PI/sH2 * pow2(alpS) / 128.0
1648     * ( 24.0 * (1.0 - 2*tHSq*uHSq/sH2) - 8.0/3.0 );
1649
1650   // Helicity-dependent factors
1651   sigma = 0.0;
1652   for (int ha=-1;ha<=1;ha += 2) {
1653     for (int hb=-1;hb<=1;hb += 2) {
1654       // Divide by 4 for helicity average      
1655       sigma += comFacHat / 4.0 
1656         * ( (1.0-ha*hb) 
1657             - 2.0 * sH*m2Sq/tHSq/uHSq 
1658             * ( 1.0 - ha*hb - sH*m2Sq/tHSq/uHSq));
1659     }    
1660   }  
1661
1662 }
1663
1664 //--------------------------------------------------------------------------
1665
1666 // Select identity, colour and anticolour.
1667
1668 void Sigma2gg2squarkantisquark::setIdColAcol() {
1669
1670   // Set flavours.
1671   setId( id1, id2, id3Sav, id4Sav);                 
1672
1673   // Set color flow (random for now)
1674   double R = rndmPtr->flat();
1675   if (R < 0.5) setColAcol( 1, 2, 2, 3, 1, 0, 0, 3);
1676   else setColAcol( 1, 2, 3, 1, 3, 0, 0, 2);
1677
1678 }
1679
1680 //==========================================================================
1681
1682 // Sigma2qg2squarkgluino
1683 // Cross section for squark-gluino production
1684
1685 //--------------------------------------------------------------------------
1686
1687 // Initialize process. 
1688   
1689 void Sigma2qg2squarkgluino::initProc() {
1690
1691   //Typecast to the correct couplings
1692   coupSUSYPtr = (CoupSUSY*) couplingsPtr;
1693
1694   // Derive name
1695   nameSave = "q g -> "+particleDataPtr->name(abs(id3Sav))+" gluino + c.c.";
1696
1697   // Final-state mass squares
1698   m2Glu     = pow2(particleDataPtr->m0(1000021));
1699   m2Sq      = pow2(particleDataPtr->m0(id3Sav));
1700
1701   // Secondary open width fraction.
1702   openFracPair = particleDataPtr->resOpenFrac(id3Sav, 1000021);
1703
1704 }
1705
1706 //--------------------------------------------------------------------------
1707
1708 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour. 
1709
1710 void Sigma2qg2squarkgluino::sigmaKin() {
1711   
1712   // Common pre-factor
1713   comFacHat = (M_PI / sH2) * pow2(alpS) * 0.5 * openFracPair;
1714
1715   // Invariants (still with Pythia 6 sign convention)
1716   double tGlu = m2Glu-tH;
1717   double uGlu = m2Glu-uH;
1718   double tSq  = m2Sq-tH;
1719   double uSq  = m2Sq-uH;
1720
1721   // Color flow A: quark color annihilates with anticolor of g
1722   sigmaA = 0.5*4./9.* tGlu/sH + (tGlu*sH+2.*m2Glu*tSq)/pow2(tGlu) -
1723     ( (sH-m2Sq+m2Glu)*(-tSq)-sH*m2Glu )/sH/(-tGlu) 
1724     + 0.5*1./2.*( tSq*(tH+2.*uH+m2Glu)-tGlu*(sH-2.*tSq) 
1725                   + (-uGlu)*(tH+m2Glu+2.*m2Sq) )/2./tGlu/uSq;
1726   // Color flow B: quark and gluon colors iterchanged
1727   sigmaB =     4./9.*(-uGlu)*(uH+m2Sq)/pow2(uSq) 
1728     + 1./18.* (sH*(uH+m2Glu) + 2.*(m2Sq-m2Glu)*uGlu)/sH/(-uSq) 
1729     + 0.5*4./9.*tGlu/sH 
1730     + 0.5*1./2.*(tSq*(tH+2.*uH+m2Glu)-tGlu*(sH-2.*tSq)
1731                  + (-uGlu)*(tH+m2Glu+2.*m2Sq))/2./tGlu/uSq;
1732
1733 }
1734
1735 double Sigma2qg2squarkgluino::sigmaHat() {
1736   
1737   // Check whether right incoming flavor
1738   int idQA = (id1 == 21) ? abs(id2) : abs(id1);
1739   int idSqSM = id3Sav%1000000;
1740   if (idQA != idSqSM) return 0.0;
1741   // NOTE: ONLY WORKS FOR SLHA1 ENUMERATION !!!
1742   // (should replace this by squark mixing matrix squares)
1743
1744   return comFacHat * (sigmaA + sigmaB);
1745
1746 }
1747
1748 //--------------------------------------------------------------------------
1749
1750 // Select identity, colour and anticolour.
1751
1752 void Sigma2qg2squarkgluino::setIdColAcol() {
1753
1754   // Check if charge conjugate final state?
1755   int idQ = (id1 == 21) ? id2 : id1;
1756   id3 = (idQ > 0) ? id3Sav : -id3Sav;
1757   id4 = 1000021;
1758   
1759   // Set flavors
1760   setId( id1, id2, id3, id4);                 
1761
1762   // Select color flow A or B (see above)
1763   double R = rndmPtr->flat()*(sigmaA+sigmaB);  
1764   if (idQ == id1) {
1765     setColAcol(1,0,2,1,3,0,2,3);
1766     if (R > sigmaA) setColAcol(1,0,2,3,2,0,1,3);
1767   } else {
1768     setColAcol(2,1,1,0,3,0,2,3);
1769     if (R > sigmaB) setColAcol(2,3,1,0,2,0,1,3);    
1770   }
1771   if (idQ < 0) swapColAcol();
1772
1773   // Use reflected kinematics if gq initial state
1774   if (id1 == 21) swapTU = true;
1775
1776 }
1777
1778 //==========================================================================
1779
1780 // Sigma2gg2gluinogluino
1781 // Cross section for gluino pair production from gg initial states
1782 // (validated against Pythia 6)
1783
1784 //--------------------------------------------------------------------------
1785
1786 // Initialize process. 
1787   
1788 void Sigma2gg2gluinogluino::initProc() {
1789
1790   //Typecast to the correct couplings
1791   coupSUSYPtr = (CoupSUSY*) couplingsPtr;
1792
1793   // Secondary open width fraction.
1794   openFracPair = particleDataPtr->resOpenFrac(1000021, 1000021);
1795   
1796
1797
1798 //--------------------------------------------------------------------------
1799
1800 // Evaluate d(sigmaHat)/d(tHat) - no incoming flavour dependence. 
1801
1802 void Sigma2gg2gluinogluino::sigmaKin() { 
1803
1804   // Modified Mandelstam variables for massive kinematics with m3 = m4.
1805   // tHG = tHat - m_gluino^2; uHG = uHat - m_gluino^2. 
1806   double s34Avg = 0.5 * (s3 + s4) - 0.25 * pow2(s3 - s4) / sH; 
1807   double tHG    = -0.5 * (sH - tH + uH);
1808   double uHG    = -0.5 * (sH + tH - uH); 
1809   double tHG2   = tHG * tHG;
1810   double uHG2   = uHG * uHG;
1811
1812   // Calculate kinematics dependence.
1813   sigTS  = (tHG * uHG - 2. * s34Avg * (tHG + 2. * s34Avg)) / tHG2
1814          + (tHG * uHG + s34Avg * (uHG - tHG)) / (sH * tHG);  
1815   sigUS  = (tHG * uHG - 2. * s34Avg * (uHG + 2. * s34Avg)) / uHG2
1816          + (tHG * uHG + s34Avg * (tHG - uHG)) / (sH * uHG);
1817   sigTU  = 2. * tHG * uHG / sH2 + s34Avg * (sH - 4. * s34Avg) 
1818          / (tHG * uHG);
1819   sigSum = sigTS + sigUS + sigTU;
1820     
1821   // Answer contains factor 1/2 from identical gluinos.
1822   sigma  = (M_PI / sH2) * pow2(alpS) * (9./4.) * 0.5 * sigSum 
1823          * openFracPair;  
1824
1825 }
1826
1827 //--------------------------------------------------------------------------
1828
1829 // Select identity, colour and anticolour.
1830
1831 void Sigma2gg2gluinogluino::setIdColAcol() {
1832
1833   // Flavours are trivial.
1834   setId( id1, id2, 1000021, 1000021);
1835
1836   // Three colour flow topologies, each with two orientations.
1837   double sigRand = sigSum * rndmPtr->flat();
1838   if (sigRand < sigTS) setColAcol( 1, 2, 2, 3, 1, 4, 4, 3);
1839   else if (sigRand < sigTS + sigUS) 
1840                        setColAcol( 1, 2, 3, 1, 3, 4, 4, 2);
1841   else                 setColAcol( 1, 2, 3, 4, 1, 4, 3, 2); 
1842   if (rndmPtr->flat() > 0.5) swapColAcol();
1843
1844 }
1845
1846 //==========================================================================
1847
1848 // Sigma2qqbar2gluinogluino
1849 // Cross section for gluino pair production from qqbar initial states
1850 // (validated against Pythia 6)
1851
1852 //--------------------------------------------------------------------------
1853
1854 // Initialize process. 
1855   
1856 void Sigma2qqbar2gluinogluino::initProc() {
1857
1858   //Typecast to the correct couplings
1859   coupSUSYPtr = (CoupSUSY*) couplingsPtr;
1860
1861   // Secondary open width fraction.
1862   openFracPair = particleDataPtr->resOpenFrac(1000021, 1000021);
1863   
1864
1865
1866 //--------------------------------------------------------------------------
1867
1868 // Begin evaluate d(sigmaHat)/d(tHat); flavour-independent part. 
1869
1870 void Sigma2qqbar2gluinogluino::sigmaKin() { 
1871
1872   // Modified Mandelstam variables for massive kinematics with m3 = m4.
1873   // tHG = tHat - m_gluino^2; uHG = uHat - m_gluino^2. 
1874   // (Note: tHG and uHG defined with opposite sign wrt Pythia 6)
1875   s34Avg = 0.5 * (s3 + s4) - 0.25 * pow2(s3 - s4) / sH; 
1876   tHG    = -0.5 * (sH - tH + uH);
1877   uHG    = -0.5 * (sH + tH - uH); 
1878   tHG2   = tHG * tHG;
1879   uHG2   = uHG * uHG;
1880
1881   // s-channel contribution.
1882   sigS   = (tHG2 + uHG2 + 2. * s34Avg * sH) / sH2; 
1883
1884 }
1885
1886 //--------------------------------------------------------------------------
1887
1888
1889 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence. 
1890
1891 double Sigma2qqbar2gluinogluino::sigmaHat() {  
1892
1893   // Squarks (L/R or 1/2) can contribute in t or u channel.
1894   // Assume identical CKM matrices in quark and squark sector. 
1895   // (Note: tHQL,R and uHQL,R defined with opposite sign wrt Pythia 6. 
1896   //  This affects the sign of the interference term below)
1897   double sQL    = pow2( particleDataPtr->m0(1000000 + abs(id1)) );
1898   double tHQL   = tHG + s34Avg - sQL; 
1899   double uHQL   = uHG + s34Avg - sQL; 
1900   double sQR    = pow2( particleDataPtr->m0(2000000 + abs(id1)) );
1901   double tHQR   = tHG + s34Avg - sQR; 
1902   double uHQR   = uHG + s34Avg - sQR; 
1903  
1904   // Calculate kinematics dependence.
1905   double sigQL  = (4./9.) * (tHG2 / pow2(tHQL) + uHG2 / pow2(uHQL)) 
1906                 + (1./9.) * s34Avg * sH / (tHQL * uHQL)
1907                 + (tHG2 + sH * s34Avg) /(sH * tHQL)   
1908                 + (uHG2 + sH * s34Avg) /(sH * uHQL);   
1909   double sigQR  = (4./9.) * (tHG2 / pow2(tHQR) + uHG2 / pow2(uHQR)) 
1910                 + (1./9.) * s34Avg * sH / (tHQR * uHQR)
1911                 + (tHG2 + sH * s34Avg) /(sH * tHQR)   
1912                 + (uHG2 + sH * s34Avg) /(sH * uHQR);
1913   double sigSum = sigS + 0.5 * (sigQL + sigQR);     
1914     
1915   // Answer contains factor 1/2 from identical gluinos.
1916   double sigma  = (M_PI / sH2) * pow2(alpS) * (8./3.) * 0.5 * sigSum 
1917                 * openFracPair;  
1918   return sigma;
1919
1920 }
1921
1922 //--------------------------------------------------------------------------
1923
1924 // Select identity, colour and anticolour.
1925
1926 void Sigma2qqbar2gluinogluino::setIdColAcol() {
1927
1928   // Flavours are trivial.
1929   setId( id1, id2, 1000021, 1000021);
1930
1931   // Two colour flow topologies. Swap if first is antiquark.
1932   if (rndmPtr->flat() < 0.5) setColAcol( 1, 0, 0, 2, 1, 3, 3, 2);
1933   else                    setColAcol( 1, 0, 0, 2, 3, 2, 1, 3); 
1934   if (id1 < 0) swapColAcol();
1935
1936 }
1937
1938 //==========================================================================
1939
1940 // Sigma1qq2antisquark
1941 // R-parity violating squark production
1942
1943 //--------------------------------------------------------------------------
1944
1945 // Initialise process
1946
1947 void Sigma1qq2antisquark::initProc(){
1948
1949   //Typecast to the correct couplings
1950   coupSUSYPtr = (CoupSUSY*) couplingsPtr;
1951
1952   //Construct name of the process from lambda'' couplings
1953
1954   nameSave = "q q' -> " + coupSUSYPtr->getName(idRes)+"* + c.c";
1955   codeSave = 2000 + 10*abs(idRes)/1000000 + abs(idRes)%10;
1956 }
1957
1958 //--------------------------------------------------------------------------
1959
1960 // Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour. 
1961
1962 void Sigma1qq2antisquark::sigmaKin() {
1963
1964   // Check if at least one RPV coupling non-zero
1965   if(!coupSUSYPtr->isUDD) {
1966     sigBW = 0.0;
1967     return;
1968   }
1969
1970   mRes = particleDataPtr->m0(abs(idRes));
1971   GammaRes = particleDataPtr->mWidth(abs(idRes));
1972   m2Res = pow2(mRes);
1973     
1974   sigBW        = sH * GammaRes/ ( pow2(sH - m2Res) + pow2(mRes * GammaRes) );
1975   sigBW       *= 2.0/3.0/mRes;
1976
1977   // Width out only includes open channels. 
1978   widthOut     = GammaRes * particleDataPtr->resOpenFrac(id3);
1979 }
1980
1981 //--------------------------------------------------------------------------
1982
1983 // Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence. 
1984
1985 double Sigma1qq2antisquark::sigmaHat() {
1986
1987   // Only allow (anti)quark-(anti)quark incoming states
1988   if (id1*id2 <= 0) return 0.0;    
1989
1990   // Generation indices
1991   int iA = (abs(id1)+1)/2;
1992   int iB = (abs(id2)+1)/2;
1993
1994   //Covert from pdg-code to ~u_i/~d_i basis
1995   bool idown = (abs(idRes)%2 == 1) ? true : false;
1996   int iC = (abs(idRes)/1000000 == 2) 
1997          ? (abs(idRes)%10+1)/2 + 3: (abs(idRes)%10+1)/2; 
1998
1999   // UDD structure
2000   if (abs(id1)%2 == 0 && abs(id2)%2 == 0) return 0.0;  
2001   if (abs(id1)%2 == 1 && abs(id2)%2 == 1 && idown) return 0.0;
2002   if ((abs(id1) + abs(id2))%2 == 1 && !idown) return 0.0;
2003
2004   double sigma = 0.0;
2005
2006   if(!idown){
2007    // d_i d_j --> ~u*_k
2008     for(int isq = 1; isq <=3; isq++){
2009       // Loop over R-type squark contributions
2010       sigma += pow2(coupSUSYPtr->rvUDD[isq][iA][iB]) 
2011         * norm(coupSUSYPtr->Rusq[iC][isq+3]);
2012     }
2013   }else{
2014     // u_i d_j --> d*_k
2015     // Pick the correct coupling for d-u in-state
2016     if(abs(id1)%2==1){
2017       iA = (abs(id2)+1)/2;
2018       iB = (abs(id1)+1)/2;
2019     }
2020     for(int isq = 1; isq <= 3; isq++){
2021       // Loop over R-type squark contributions
2022       sigma += pow2(coupSUSYPtr->rvUDD[iA][iB][isq]) 
2023         * norm(coupSUSYPtr->Rdsq[iC][isq+3]);
2024     }
2025   }
2026
2027   sigma *= sigBW;
2028   return sigma;    
2029
2030 }
2031
2032 //--------------------------------------------------------------------------
2033
2034 // Select identity, colour and anticolour.
2035
2036 void Sigma1qq2antisquark::setIdColAcol() {
2037
2038   // Set flavours.
2039   if(id1 < 0 && id2 < 0 ) setId( id1, id2, idRes);
2040   else setId( id1, id2, -idRes);
2041
2042   // Colour flow topologies. Swap when antiquarks.
2043   if (abs(id1) < 9) setColAcol( 1, 0, 2, 0, 0, 3);
2044   else              setColAcol( 0, 0, 0, 0, 0, 0);
2045   if (id1 < 0) swapColAcol();
2046
2047 }
2048
2049
2050 //==========================================================================
2051
2052 } // end namespace Pythia8
2053