using option '-treename HLTesdTree' for EsdCollector, adding default parameter for...
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8130 / src / SigmaSUSY.cxx
CommitLineData
5ad4eb21 1// SigmaSUSY.cc is a part of the PYTHIA event generator.
2// Copyright (C) 2008 Torbjorn Sjostrand.
3// PYTHIA is licenced under the GNU GPL version 2, see COPYING for details.
4// Please respect the MCnet Guidelines, see GUIDELINES for details.
5
6// Function definitions (not found in the header) for the
7// supersymmetry simulation classes.
8
9#include "SigmaSUSY.h"
10
11namespace Pythia8 {
12
13//**************************************************************************
14
15// Sigma2qqbar2gauginogaugino "mother" class.
16// Cross section for gaugino pair production: neutralino pair,
17// neutralino-chargino, and chargino pair production all inherit from this.
18
19//*********
20
21// Initialize process.
22
23void Sigma2qqbar2gauginogaugino::initProc() {
24
25 // Construct name of process.
26 nameSave = "q qbar -> " + ParticleDataTable::name(id3) + " "
27 + ParticleDataTable::name(id4);
28
29 // Count number of final-state charged particles (default = 0)
30 // (useful since charginos inherit from this class.)
31 nCharged=0;
32 if (abs(id3) == 1000024 || abs(id3) == 1000037) nCharged++;
33 if (abs(id4) == 1000024 || abs(id4) == 1000037) nCharged++;
34
35 // Set up couplings
36 mZpole = ParticleDataTable::m0(23);
37 wZpole = ParticleDataTable::mWidth(23);
38 double mWpole = ParticleDataTable::m0(24);
39
40 // Running masses and weak mixing angle
41 // (default to pole values if no running available)
42 double mW = mWpole;
43 double mZ = mZpole;
44 sin2W = 1.0 - pow(mW/mZ,2);
45 if (slhaPtr->gauge.exists(1) && slhaPtr->gauge.exists(2)
46 && slhaPtr->hmix.exists(3)) {
47 double gp=slhaPtr->gauge(1);
48 double g =slhaPtr->gauge(2);
49 double v =slhaPtr->hmix(3);
50 mW = g * v / 2.0;
51 mZ = sqrt(pow(gp,2)+pow(g,2)) * v / 2.0;
52 sin2W = pow2(gp)/(pow2(g)+pow2(gp));
53 }
54
55 // Shorthand for SUSY couplings
56 // By default, use the running one in HMIX
57 // If not found, use the MZ one in MINPAR
58 double tanb = slhaPtr->hmix.exists(2) ?
59 slhaPtr->hmix(2) : slhaPtr->minpar(3);
60 double sinW=sqrt(sin2W);
61 double cosW=sqrt(1.0-sin2W);
62 double cosb = sqrt( 1.0 / (1.0 + tanb*tanb) );
63 double sinb = sqrt(max(0.0,1.0-cosb*cosb));
64 SusyLesHouches::matrixblock<6> Ru(slhaPtr->usqmix);
65 SusyLesHouches::matrixblock<6> Rd(slhaPtr->dsqmix);
66 SusyLesHouches::matrixblock<6> imRu(slhaPtr->imusqmix);
67 SusyLesHouches::matrixblock<6> imRd(slhaPtr->imusqmix);
68
69 // Local complex copies of neutralino mixing matrix entries.
70 // (incl NMSSM for future use)
71 complex ni1,nj1,ni2,nj2,ni3,nj3,ni4,nj4,ni5,nj5;
72 if (slhaPtr->modsel(3) != 1) {
73 ni1=complex( slhaPtr->nmix(id3chi,1), slhaPtr->imnmix(id3chi,1) );
74 nj1=complex( slhaPtr->nmix(id4chi,1), slhaPtr->imnmix(id4chi,1) );
75 ni2=complex( slhaPtr->nmix(id3chi,2), slhaPtr->imnmix(id3chi,2) );
76 nj2=complex( slhaPtr->nmix(id4chi,2), slhaPtr->imnmix(id4chi,2) );
77 ni3=complex( slhaPtr->nmix(id3chi,3), slhaPtr->imnmix(id3chi,3) );
78 nj3=complex( slhaPtr->nmix(id4chi,3), slhaPtr->imnmix(id4chi,3) );
79 ni4=complex( slhaPtr->nmix(id3chi,4), slhaPtr->imnmix(id3chi,4) );
80 nj4=complex( slhaPtr->nmix(id4chi,4), slhaPtr->imnmix(id4chi,4) );
81 ni5=complex( 0.,0.);
82 nj5=complex( 0.,0.);
83 } else {
84 ni1=complex( slhaPtr->nmnmix(id3chi,1), slhaPtr->imnmnmix(id3chi,1) );
85 nj1=complex( slhaPtr->nmnmix(id4chi,1), slhaPtr->imnmnmix(id4chi,1) );
86 ni2=complex( slhaPtr->nmnmix(id3chi,2), slhaPtr->imnmnmix(id3chi,2) );
87 nj2=complex( slhaPtr->nmnmix(id4chi,2), slhaPtr->imnmnmix(id4chi,2) );
88 ni3=complex( slhaPtr->nmnmix(id3chi,3), slhaPtr->imnmnmix(id3chi,3) );
89 nj3=complex( slhaPtr->nmnmix(id4chi,3), slhaPtr->imnmnmix(id4chi,3) );
90 ni4=complex( slhaPtr->nmnmix(id3chi,4), slhaPtr->imnmnmix(id3chi,4) );
91 nj4=complex( slhaPtr->nmnmix(id4chi,4), slhaPtr->imnmnmix(id4chi,4) );
92 ni5=complex( slhaPtr->nmnmix(id3chi,5), slhaPtr->imnmnmix(id3chi,5) );
93 nj5=complex( slhaPtr->nmnmix(id4chi,5), slhaPtr->imnmnmix(id4chi,5) );
94 }
95
96 // Change to positive mass convention.
97 complex iRot( 0., 1.);
98 if (slhaPtr->mass(abs(id3)) < 0.) {
99 ni1 *= iRot;
100 ni2 *= iRot;
101 ni3 *= iRot;
102 ni4 *= iRot;
103 ni5 *= iRot;
104 };
105 if (slhaPtr->mass(abs(id4)) < 0.) {
106 nj1 *= iRot;
107 nj2 *= iRot;
108 nj3 *= iRot;
109 nj4 *= iRot;
110 nj5 *= iRot;
111 };
112
113 // Local copies of Chargino mixing
114 complex ui1,ui2,vi1,vi2,uj1,uj2,vj1,vj2;
115 if (id3chi <= 2) {
116 ui1=complex( slhaPtr->umix(id3chi,1), slhaPtr->imumix(id3chi,1) );
117 ui2=complex( slhaPtr->umix(id3chi,2), slhaPtr->imumix(id3chi,2) );
118 vi1=complex( slhaPtr->vmix(id3chi,1), slhaPtr->imvmix(id3chi,1) );
119 vi2=complex( slhaPtr->vmix(id3chi,2), slhaPtr->imvmix(id3chi,2) );
120 }
121 if (id4chi <= 2) {
122 uj1=complex( slhaPtr->umix(id4chi,1), slhaPtr->imumix(id4chi,1) );
123 uj2=complex( slhaPtr->umix(id4chi,2), slhaPtr->imumix(id4chi,2) );
124 vj1=complex( slhaPtr->vmix(id4chi,1), slhaPtr->imvmix(id4chi,1) );
125 vj2=complex( slhaPtr->vmix(id4chi,2), slhaPtr->imvmix(id4chi,2) );
126 }
127
128 // Z chi_i chi_j
129 OLpp = -0.5 * ni3 * conj(nj3) + 0.5 * ni4 * conj(nj4);
130 ORpp = 0.5 * conj(ni3) * nj3 - 0.5 * conj(ni4) * nj4;
131
132 // Z cha_i cha_j
133 OLp = -vi1*conj(vj1) - 0.5*vi2*conj(vj2)
134 + ( (id3chi == id4chi) ? sin2W : 0.0);
135 ORp = -conj(ui1)*uj1 - 0.5*conj(ui2)*uj2
136 + ( (id3chi == id4chi) ? sin2W : 0.0);
137
138 // W chi_i cha_j
139 OL = -1.0/sqrt(2.0)*ni4*conj(vj2)+ni2*conj(vj1);
140 OR = 1.0/sqrt(2.0)*conj(ni3)*uj2+conj(ni2)*uj1;
141
142 // Z q_{idq} q_{idq} (def with extra factor 2 compared to [Okun])
143 for (int idq = 1; idq <= 5; ++idq) {
144 // No FCNC in Zqq, so here sum over diagonal only
145 // No t in beam, so only sum up to 5
146 LqqZ[idq] = CoupEW::lf(idq);
147 RqqZ[idq] = CoupEW::rf(idq);
148 }
149
150 // ~chi^0_i ~q_jsq idq
151 for (int jsq = 1; jsq<=6; jsq++) {
152 // Sum jsq over all squarks
153 for (int idq = 1; idq<=5; idq++) {
154 // No t in beam, so only sum iq over 5 lightest quarks.
155
156 // quark index
157 int k=(idq+1)/2;
158
159 // Set quark mass for use in couplings
160 // Initial guess 0,0,0,mc,mb, with the latter from the PDT
161 double mq = ParticleDataTable::m0(idq);
162 if (idq <= 3) mq=0.0;
163
164 // Treat u and d quarks separately
165 if (idq % 2 == 1) {
166
167 // idq = d quark
168 double eq = -1.0/3.0;
169 double T3q = -0.5;
170
171 // Compute running mass from Yukawas and vevs if possible.
172 if (slhaPtr->yd.exists() && slhaPtr->hmix.exists(3)) {
173 double ykk=slhaPtr->yd(k,k);
174 double v1=slhaPtr->hmix(3)/sqrt(1+pow(tanb,2));
175 if (ykk != 0.0) mq = ykk * v1 / sqrt(2.0) ;
176 // cout <<scientific<<" q "<<idq << " (k="<<k<<") Y = "
177 // <<ykk<<" v = "<<v1<<" m = "<<mq<<endl;
178 }
179
180 // Shorthand for squark mixing matrices
181 complex RdjkL(Rd(jsq,k),imRd(jsq,k));
182 complex RdjkR(Rd(jsq,k+3),imRd(jsq,k+3));
183
184 // L(~d_jsq,d_k,~chi^0_i)
185 LsqXi[jsq][idq]=((eq-T3q)*sinW*ni1+T3q*cosW*ni2) * conj(RdjkL)
186 + mq*cosW*ni3*conj(RdjkR)/2.0/mW/cosb ;
187
188 // L(~d_jsq,d_k,~chi^0_j)
189 LsqXj[jsq][idq]=((eq-T3q)*sinW*nj1+T3q*cosW*nj2) * conj(RdjkL)
190 + mq*cosW*nj3*conj(RdjkR)/2.0/mW/cosb ;
191
192 // R(~d_jsq,d_k,~chi^0_i)
193 RsqXi[jsq][idq]=eq*sinW*ni1*RdjkR - mq*cosW*ni3*RdjkL/2.0/mW/cosb ;
194 RsqXi[jsq][idq]=-conj(RsqXi[jsq][idq]) ;
195
196 // R(~d_jsq,d_k,~chi^0_j)
197 RsqXj[jsq][idq]=eq*sinW*nj1*RdjkR - mq*cosW*nj3*RdjkL/2.0/mW/cosb ;
198 RsqXj[jsq][idq]=-conj(RsqXj[jsq][idq]) ;
199
200 // Chargino couplings (initialize)
201 LsqCi[jsq][idq]=0.0;
202 RsqCi[jsq][idq]=0.0;
203 LsqCj[jsq][idq]=0.0;
204 RsqCj[jsq][idq]=0.0;
205
206 // Sum over flavours (off-diagonal mixings)
207 // (note l <-> k with respect to [Klasen et al])
208 for (int l=1; l<=3; l++) {
209
210 // Set quark masses for use in couplings
211 // Initial guess 0,0,0,mc,mb, with the latter from the PDT
212 double mu = ParticleDataTable::m0(2*l);
213 if (2*l <= 3) mu=0.0;
214
215 // Compute running u mass from Yukawas and vevs if possible.
216 if (slhaPtr->yu.exists() && slhaPtr->hmix.exists(3)) {
217 double yll=slhaPtr->yu(l,l);
218 double v2=slhaPtr->hmix(3)/sqrt(1.0+1.0/pow(tanb,2));
219 if (yll != 0.0) mu = yll * v2 / sqrt(2.0) ;
220 }
221
222 // Shorthand for u squark mixing
223 complex RujlL(Ru(jsq,l),imRu(jsq,l));
224 complex RujlR(Ru(jsq,l+3),imRu(jsq,l+3));
225
226 // CKM matrix (use Pythia-8 one if no SLHA)
227 complex Vlk=VCKM::Vid(l,k);
228 if (slhaPtr->vckm.exists())
229 Vlk=complex(slhaPtr->vckm(l,k),slhaPtr->imvckm(l,k));
230
231 // L(~u_jsq,d_k,~chi^+/-_i)
232 LsqCi[jsq][idq] += Vlk*(conj(vi1)*RujlL
233 -mu*conj(vi2)*RujlR/sqrt(2.0)/mW/sinb );
234
235 // L(~u_jsq,d_k,~chi^+/-_j)
236 LsqCj[jsq][idq] += Vlk*(conj(vj1)*RujlL
237 -mu*conj(vj2)*RujlR/sqrt(2.0)/mW/sinb );
238
239 // R(~u_jsq,d_k,~chi^+/-_i)
240 RsqCi[jsq][idq] += - Vlk*mq*ui2*RujlL/sqrt(2.0)/mW/cosb ;
241
242 // R(~u_jsq,d_k,~chi^+/-_j)
243 RsqCj[jsq][idq] += - Vlk*mq*uj2*RujlL/sqrt(2.0)/mW/cosb ;
244
245 }
246
247 } else {
248
249 // idq = u quark
250 double eq = 2.0/3.0;
251 double T3q = 0.5;
252
253 // Compute mass from Yukawas and vevs if possible.
254 if (slhaPtr->yu.exists() && slhaPtr->hmix.exists(3)) {
255 double ykk=slhaPtr->yu(k,k);
256 double v2=slhaPtr->hmix(3)/sqrt(1.0+1.0/pow(tanb,2));
257 if (ykk != 0.0) mq = ykk * v2 / sqrt(2.0) ;
258 // cout <<scientific<<" q "<<idq << " (k="<<k<<") Y = "
259 // <<ykk<<" v = "<<v2<<" m = "<<mq<<endl;
260 }
261
262 // Shorthand for squark mixing matrices
263 complex RujkL(Ru(jsq,k),imRu(jsq,k));
264 complex RujkR(Ru(jsq,k+3),imRu(jsq,k+3));
265
266 // L(~u_jsq,u_k,~chi^0_i)
267 LsqXi[jsq][idq]=((eq-T3q)*sinW*ni1+T3q*cosW*ni2) * conj(RujkL)
268 + mq*cosW*ni4*conj(RujkR)/2.0/mW/sinb ;
269
270 // L(~u_jsq,u_k,~chi^0_j)
271 LsqXj[jsq][idq]=((eq-T3q)*sinW*nj1+T3q*cosW*nj2) * conj(RujkL)
272 + mq*cosW*nj4*conj(RujkR)/2.0/mW/sinb ;
273
274 // R(~u_jsq,u_k,~chi^0_i)
275 RsqXi[jsq][idq]=eq*sinW*ni1*RujkR - mq*cosW*ni4*RujkL/2.0/mW/sinb ;
276 RsqXi[jsq][idq]=-conj(RsqXi[jsq][idq]);
277
278 // R(~u_jsq,u_k,~chi^0_i)
279 RsqXj[jsq][idq]=eq*sinW*nj1*RujkR - mq*cosW*nj4*RujkL/2.0/mW/sinb ;
280 RsqXj[jsq][idq]=-conj(RsqXj[jsq][idq]);
281
282 // Chargino couplings (initialize)
283 LsqCi[jsq][idq]=0.0;
284 RsqCi[jsq][idq]=0.0;
285 LsqCj[jsq][idq]=0.0;
286 RsqCj[jsq][idq]=0.0;
287
288 // Sum over flavours (off-diagonal mixings)
289 // (note l <-> k with respect to [Klasen et al])
290 for (int l=1; l<=3; l++) {
291
292 // Set quark masses for use in couplings
293 // Initial guess 0,0,0,mc,mb, with the latter from the PDT
294 double md = ParticleDataTable::m0(2*l-1);
295 if (2*l-1 <= 3) md=0.0;
296
297 // Compute running d mass from Yukawas and vevs if possible.
298 if (slhaPtr->yd.exists() && slhaPtr->hmix.exists(3)) {
299 double yll=slhaPtr->yd(l,l);
300 double v1=slhaPtr->hmix(3)/sqrt(1+pow(tanb,2));
301 if (yll != 0.0) md = yll * v1 / sqrt(2.0) ;
302 }
303
304 // Shorthand for d squark mixing
305 complex RdjlL(Rd(jsq,l),imRd(jsq,l));
306 complex RdjlR(Rd(jsq,l+3),imRd(jsq,l+3));
307
308 // CKM matrix (use Pythia-8 one if no SLHA)
309 complex Vkl=VCKM::Vid(k,l);
310 if (slhaPtr->vckm.exists())
311 Vkl=complex(slhaPtr->vckm(k,l),slhaPtr->imvckm(k,l));
312
313 // L(~d_jsq,u_k,~chi^+/-_i)
314 LsqCi[jsq][idq] += Vkl*(ui1*conj(RdjlL)
315 -md*ui2*conj(RdjlR)/sqrt(2.0)/mW/cosb) ;
316
317 // L(~d_jsq,u_k,~chi^+/-_j)
318 LsqCj[jsq][idq] += Vkl*(uj1*conj(RdjlL)
319 -md*uj2*conj(RdjlR)/sqrt(2.0)/mW/cosb) ;
320
321 // R(~d_jsq,u_k,~chi^+/-_i)
322 RsqCi[jsq][idq] += - Vkl*mq*conj(vi2*RdjlL)/sqrt(2.0)/mW/sinb ;
323
324 // R(~d_jsq,u_k,~chi^+/-_j)
325 RsqCj[jsq][idq] += - Vkl*mq*conj(vi2*RdjlL)/sqrt(2.0)/mW/sinb ;
326
327 }
328
329 }
330 }
331 }
332}
333
334//*********
335
336// Evaluate d(sigmaHat)/d(tHat), part independent of incoming flavour.
337
338void Sigma2qqbar2gauginogaugino::sigmaKin() {
339
340 // Common flavour-independent factor.
341 sigma0 = (M_PI / sH2) * pow2(alpEM) ;
342
343 // Factor 1/2 for identical final particles.
344 if (id3 == id4 && nCharged == 0) sigma0 *= 0.5;
345
346 // Auxiliary factors for use below
347 ui = uH - s3;
348 uj = uH - s4;
349 ti = tH - s3;
350 tj = tH - s4;
351 sz = sH - pow2(mZpole);
352 d = pow2(sz) + pow2(mZpole * wZpole);
353 propZ = complex( sz / d, mZpole * wZpole / d);
354
355}
356
357//*********
358
359// Evaluate d(sigmaHat)/d(tHat), including incoming flavour dependence.
360
361// neutralino pairs:
362
363double Sigma2qqbar2gauginogaugino::sigmaHat() {
364
365 // Only allow quark-antiquark incoming states
366 if (id1*id2 >= 0) {
367 return 0.0;
368 }
369
370 // Only allow incoming states with sum(charge) = 0
371 if ((id1+id2) % 2 != 0) {
372 return 0.0;
373 }
374
375 // Flavour-dependent kinematics-dependent couplings.
376 int idAbs1 = abs(id1);
377 int idAbs2 = abs(id2);
378 complex QuLL = (id1 == -id2) ? LqqZ[idAbs1] * OLpp/2.0 * propZ : 0.0;
379 complex QtLL = (id1 == -id2) ? LqqZ[idAbs1] * ORpp/2.0 * propZ : 0.0;
380 complex QuRR = (id1 == -id2) ? RqqZ[idAbs1] * ORpp/2.0 * propZ : 0.0;
381 complex QtRR = (id1 == -id2) ? RqqZ[idAbs1] * OLpp/2.0 * propZ : 0.0;
382 complex QuLR = 0.0;
383 complex QtLR = 0.0;
384 complex QuRL = 0.0;
385 complex QtRL = 0.0;
386
387 // Add t-channel squark flavour sums to QmXY couplings
388 for (int jsq=1; jsq<=6; jsq++) {
389 int idsq=((jsq+2)/3)*1000000 + 2*((jsq-1) % 3) + (idAbs1+1) % 2 + 1;
390 double msq2=pow(ParticleDataTable::m0(idsq),2);
391 double usq = uH - msq2;
392 double tsq = tH - msq2;
393 QuLL += LsqXi[jsq][idAbs2]*conj(LsqXj[jsq][idAbs1])/usq;
394 QtLL -= conj(LsqXi[jsq][idAbs1])*LsqXj[jsq][idAbs2]/tsq;
395 QuRR += RsqXi[jsq][idAbs2]*conj(RsqXj[jsq][idAbs1])/usq;
396 QtRR -= conj(RsqXi[jsq][idAbs1])*RsqXj[jsq][idAbs2]/tsq;
397 QuLR += RsqXi[jsq][idAbs2]*conj(LsqXj[jsq][idAbs1])/usq;
398 QtLR += conj(LsqXi[jsq][idAbs1])*RsqXj[jsq][idAbs2]/tsq;
399 QuRL += LsqXi[jsq][idAbs2]*conj(RsqXj[jsq][idAbs1])/usq;
400 QtRL += conj(RsqXi[jsq][idAbs1])*LsqXj[jsq][idAbs2]/tsq;
401 }
402
403 // Multiply by overall normalization
404 QuLL *= 1.0/(sin2W * (1 - sin2W));
405 QtLL *= 1.0/(sin2W * (1 - sin2W));
406 QuRR *= 1.0/(sin2W * (1 - sin2W));
407 QtRR *= 1.0/(sin2W * (1 - sin2W));
408 QuLR *= 1.0/(sin2W * (1 - sin2W));
409 QtLR *= 1.0/(sin2W * (1 - sin2W));
410 QuRL *= 1.0/(sin2W * (1 - sin2W));
411 QtRL *= 1.0/(sin2W * (1 - sin2W));
412
413 // Compute matrix element weight
414 double weight = 0;
415 // Average over separate helicity contributions
416 // LL (ha = -1, hb = +1) (divided by 4 for average)
417 weight += norm(QuLL) * ui * uj + norm(QtLL) * ti * tj
418 + 2 * real(conj(QuLL) * QtLL) * m3 * m4 * sH;
419 // RR (ha = 1, hb = -1) (divided by 4 for average)
420 weight += norm(QtRR) * ti * tj + norm(QuRR) * ui * uj
421 + 2 * real(conj(QuRR) * QtRR) * m3 * m4 * sH;
422 // RL (ha = 1, hb = 1) (divided by 4 for average)
423 weight += norm(QuRL) * ui * uj + norm(QtRL) * ti * tj
424 - real(conj(QuRL) * QtRL) * (uH * tH - s3 * s4);
425 // LR (ha = -1, hb = -1) (divided by 4 for average)
426 weight += norm(QuLR) * ui * uj + norm(QtLR) * ti * tj
427 - real(conj(QuLR) * QtLR) * (uH * tH - s3 * s4);
428
429 // Cross section, including colour factor.
430 double sigma = sigma0 * weight;
431 if (idAbs1 < 9) sigma /= 3.;
432
433 // Answer.
434 return sigma;
435
436}
437
438// neutralino-chargino
439
440double Sigma2qqbar2chi0char::sigmaHat() {
441
442 // This cross section not fully debugged yet, skip.
443 return 0.0;
444
445 // Only allow quark-antiquark incoming states
446 if (id1*id2 >= 0) {
447 return 0.0;
448 }
449
450 // Only allow incoming states with sum(charge) = +/- 1
451 if ((id1+id2) % 2 == 0) {
452 return 0.0;
453 }
454
455 // Flavour-dependent kinematics-dependent couplings.
456 int idAbs1 = abs(id1);
457 int idAbs2 = abs(id2);
458 complex QuLL = (id1 == -id2) ? LqqZ[idAbs1] * OLpp/2.0 * propZ : 0.0;
459 complex QtLL = (id1 == -id2) ? LqqZ[idAbs1] * ORpp/2.0 * propZ : 0.0;
460 complex QuRR = (id1 == -id2) ? RqqZ[idAbs1] * ORpp/2.0 * propZ : 0.0;
461 complex QtRR = (id1 == -id2) ? RqqZ[idAbs1] * OLpp/2.0 * propZ : 0.0;
462 complex QuLR = 0.0;
463 complex QtLR = 0.0;
464 complex QuRL = 0.0;
465 complex QtRL = 0.0;
466
467 // Add t-channel squark flavour sums to QmXY couplings
468 for (int jsq=1; jsq<=6; jsq++) {
469 int idsq=((jsq+2)/3)*1000000 + 2*((jsq-1) % 3) + (idAbs1+1) % 2 + 1;
470 double msq2=pow(ParticleDataTable::m0(idsq),2);
471 double usq = uH - msq2;
472 double tsq = tH - msq2;
473 QuLL += LsqXi[jsq][idAbs2]*conj(LsqXj[jsq][idAbs1])/usq;
474 QtLL -= conj(LsqXi[jsq][idAbs1])*LsqXj[jsq][idAbs2]/tsq;
475 QuRR += RsqXi[jsq][idAbs2]*conj(RsqXj[jsq][idAbs1])/usq;
476 QtRR -= conj(RsqXi[jsq][idAbs1])*RsqXj[jsq][idAbs2]/tsq;
477 QuLR += RsqXi[jsq][idAbs2]*conj(LsqXj[jsq][idAbs1])/usq;
478 QtLR += conj(LsqXi[jsq][idAbs1])*RsqXj[jsq][idAbs2]/tsq;
479 QuRL += LsqXi[jsq][idAbs2]*conj(RsqXj[jsq][idAbs1])/usq;
480 QtRL += conj(RsqXi[jsq][idAbs1])*LsqXj[jsq][idAbs2]/tsq;
481 }
482
483 // Compute matrix element weight
484 double weight = 0;
485 // Average over separate helicity contributions
486 // LL (ha = -1, hb = +1) (divided by 4 for average)
487 weight += norm(QuLL) * ui * uj + norm(QtLL) * ti * tj
488 + 2 * real(conj(QuLL) * QtLL) * m3 * m4 * sH;
489 // RR (ha = 1, hb = -1) (divided by 4 for average)
490 weight += norm(QtRR) * ti * tj + norm(QuRR) * ui * uj
491 + 2 * real(conj(QuRR) * QtRR) * m3 * m4 * sH;
492 // RL (ha = 1, hb = 1) (divided by 4 for average)
493 weight += norm(QuRL) * ui * uj + norm(QtRL) * ti * tj
494 - real(conj(QuRL) * QtRL) * (uH * tH - s3 * s4);
495 // LR (ha = -1, hb = -1) (divided by 4 for average)
496 weight += norm(QuLR) * ui * uj + norm(QtLR) * ti * tj
497 - real(conj(QuLR) * QtLR) * (uH * tH - s3 * s4);
498
499 // Cross section, including colour factor.
500 double sigma = sigma0 * weight;
501 if (idAbs1 < 9) sigma /= 3.;
502
503 // Answer.
504 return sigma;
505
506}
507
508// chargino-chargino
509
510double Sigma2qqbar2charchar::sigmaHat() {
511
512 // This cross section not fully debugged yet, skip.
513 return 0.0;
514
515 // Only allow quark-antiquark incoming states
516 if (id1*id2 >= 0) {
517 return 0.0;
518 }
519
520 // Only allow incoming states with sum(charge) = 0
521 if ((id1+id2) % 2 != 0) {
522 return 0.0;
523 }
524
525 // Flavour-dependent kinematics-dependent couplings.
526 int idAbs1 = abs(id1);
527 int idAbs2 = abs(id2);
528
529 // Add s-channel Z
530 complex QuLL = (id1 == -id2) ? - LqqZ[idAbs1] * conj(OLp)/2.0 * propZ
531 / (sin2W * (1 - sin2W)) : 0.0;
532 complex QtLL = (id1 == -id2) ? - LqqZ[idAbs1] * conj(ORp)/2.0 * propZ
533 / (sin2W * (1 - sin2W)) : 0.0;
534 complex QuRR = (id1 == -id2) ? - RqqZ[idAbs1] * conj(ORp)/2.0 * propZ
535 / (sin2W * (1 - sin2W)) : 0.0;
536 complex QtRR = (id1 == -id2) ? - RqqZ[idAbs1] * conj(OLp)/2.0 * propZ
537 / (sin2W * (1 - sin2W)) : 0.0;
538 complex QuLR = 0.0;
539 complex QtLR = 0.0;
540 complex QuRL = 0.0;
541 complex QtRL = 0.0;
542
543 // Add s-channel photon
544 if (idAbs1 == idAbs2 && abs(id3) == abs(id4)) {
545 double e = -1.0/3.0;
546 if (idAbs1 % 2 == 0) e = 2.0/3.0;
547 if (idAbs1 > 10) e = -1.0;
548 QuLL += e/sH;
549 QtLL += e/sH;
550 QuRR += e/sH;
551 QtRR += e/sH;
552 }
553
554 // Add t-channel squark flavour sums to QmXY couplings
555 for (int jsq=1; jsq<=6; jsq++) {
556 int idsq=((jsq+2)/3)*1000000 + 2*((jsq-1) % 3) + (idAbs1 % 2) + 1;
557 double msq2=pow(ParticleDataTable::m0(idsq),2);
558 double usq = uH - msq2;
559 double tsq = tH - msq2;
560
561 // up-type quarks get u-channel d-squark contributions
562 if (idAbs1 % 2 == 0 ) {
563
564 // Flip 1<->2 and i<->j if id1 is antiquark
565 if (id1 > 0) {
566 QuLL += LsqCi[jsq][idAbs2]*conj(LsqCj[jsq][idAbs1])/usq/2.0/sin2W;
567 QuRR += RsqCi[jsq][idAbs2]*conj(RsqCj[jsq][idAbs1])/usq/2.0/sin2W;
568 QuLR += RsqCi[jsq][idAbs2]*conj(LsqCj[jsq][idAbs1])/usq/2.0/sin2W;
569 QuRL += LsqCi[jsq][idAbs2]*conj(RsqCj[jsq][idAbs1])/usq/2.0/sin2W;
570 } else {
571 QuLL += LsqCj[jsq][idAbs1]*conj(LsqCi[jsq][idAbs2])/usq/2.0/sin2W;
572 QuRR += RsqCj[jsq][idAbs1]*conj(RsqCi[jsq][idAbs2])/usq/2.0/sin2W;
573 QuLR += RsqCj[jsq][idAbs1]*conj(LsqCi[jsq][idAbs2])/usq/2.0/sin2W;
574 QuRL += LsqCj[jsq][idAbs1]*conj(RsqCi[jsq][idAbs2])/usq/2.0/sin2W;
575 }
576 }
577
578 // down-type quarks get t-channel u-squark contributions
579 else if (idAbs1 % 2 == 1) {
580 // Flip 1<->2 and i<->j if id1 is antiquark
581 if (id1 > 0) {
582 QtLL -= LsqCi[jsq][idAbs1]*conj(LsqCj[jsq][idAbs2])/tsq/2.0/sin2W;
583 QtRR -= RsqCi[jsq][idAbs1]*conj(RsqCj[jsq][idAbs2])/tsq/2.0/sin2W;
584 QtLR += LsqCi[jsq][idAbs1]*conj(RsqCj[jsq][idAbs2])/tsq/2.0/sin2W;
585 QtRL += RsqCi[jsq][idAbs1]*conj(LsqCj[jsq][idAbs2])/tsq/2.0/sin2W;
586 } else {
587 QtLL -= LsqCj[jsq][idAbs2]*conj(LsqCi[jsq][idAbs1])/tsq/2.0/sin2W;
588 QtRR -= RsqCj[jsq][idAbs2]*conj(RsqCi[jsq][idAbs1])/tsq/2.0/sin2W;
589 QtLR += LsqCj[jsq][idAbs2]*conj(RsqCi[jsq][idAbs1])/tsq/2.0/sin2W;
590 QtRL += RsqCj[jsq][idAbs2]*conj(LsqCi[jsq][idAbs1])/tsq/2.0/sin2W;
591 }
592
593 }
594 }
595
596 // Compute matrix element weight
597 double weight = 0;
598 // Average over separate helicity contributions
599 // LL (ha = -1, hb = +1) (divided by 4 for average)
600 weight += norm(QuLL) * ui * uj + norm(QtLL) * ti * tj
601 + 2 * real(conj(QuLL) * QtLL) * m3 * m4 * sH;
602 // RR (ha = 1, hb = -1) (divided by 4 for average)
603 weight += norm(QtRR) * ti * tj + norm(QuRR) * ui * uj
604 + 2 * real(conj(QuRR) * QtRR) * m3 * m4 * sH;
605 // RL (ha = 1, hb = 1) (divided by 4 for average)
606 weight += norm(QuRL) * ui * uj + norm(QtRL) * ti * tj
607 - real(conj(QuRL) * QtRL) * (uH * tH - s3 * s4);
608 // LR (ha = -1, hb = -1) (divided by 4 for average)
609 weight += norm(QuLR) * ui * uj + norm(QtLR) * ti * tj
610 - real(conj(QuLR) * QtLR) * (uH * tH - s3 * s4);
611
612 // Cross section, including colour factor.
613 double sigma = sigma0 * weight;
614 if (idAbs1 < 9) sigma /= 3.;
615
616 // Answer.
617 return sigma;
618
619}
620
621//*********
622
623// Select identity, colour and anticolour.
624
625void Sigma2qqbar2gauginogaugino::setIdColAcol() {
626
627 // Set flavours.
628
629 if (nCharged == 0) {
630 // Neutralino-Neutralino
631 setId( id1, id2, id3, id4);
632
633 } else if (nCharged == 1) {
634 // Neutralino-Chargino: set charge according to in-states
635 if ( abs(id1%1) == 1 && id1 > 0 ) {
636 setId( id1, id2, id3, -id4);
637 } else {
638 setId( id1, id2, id3, id4);
639 }
640
641 } else if (nCharged == 2) {
642 // Chargino-Chargino: charge already included in id3, id4
643 setId( id1, id2, id3, id4);
644 }
645
646 // Colour flow topologies. Swap when antiquarks.
647 if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
648 else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
649 if (id1 < 0) swapColAcol();
650
651}
652
653//**************************************************************************
654
655} // end namespace Pythia8
656