StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
SigmaLeftRightSym.cc
1 // SigmaLeftRightSym.cc is a part of the PYTHIA event generator.
2 // Copyright (C) 2012 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 // left-right-symmetry simulation classes.
8 
9 #include "SigmaLeftRightSym.h"
10 
11 namespace Pythia8 {
12 
13 //==========================================================================
14 
15 // Sigma1ffbar2ZRight class.
16 // Cross section for f fbar -> Z_R^0 (righthanded gauge boson).
17 
18 //--------------------------------------------------------------------------
19 
20 // Initialize process.
21 
22 void Sigma1ffbar2ZRight::initProc() {
23 
24  // Store Z_R mass and width for propagator.
25  idZR = 9900023;
26  mRes = particleDataPtr->m0(idZR);
27  GammaRes = particleDataPtr->mWidth(idZR);
28  m2Res = mRes*mRes;
29  GamMRat = GammaRes / mRes;
30  sin2tW = couplingsPtr->sin2thetaW();
31 
32  // Set pointer to particle properties and decay table.
33  ZRPtr = particleDataPtr->particleDataEntryPtr(idZR);
34 
35 }
36 
37 //--------------------------------------------------------------------------
38 
39 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
40 
41 void Sigma1ffbar2ZRight::sigmaKin() {
42 
43  // Set up Breit-Wigner. Width out only includes open channels.
44  double sigBW = 12. * M_PI/ ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
45  double widthOut = ZRPtr->resWidthOpen(idZR, mH);
46 
47  // Prefactor for incoming widths. Combine. Done.
48  double preFac = alpEM * mH / ( 48. * sin2tW * (1. - sin2tW)
49  * (1. - 2. * sin2tW) );
50  sigma0 = preFac * sigBW * widthOut;
51 
52 }
53 
54 //--------------------------------------------------------------------------
55 
56 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
57 
58 double Sigma1ffbar2ZRight::sigmaHat() {
59 
60  // Vector and axial couplings of incoming fermion pair.
61  int idAbs = abs(id1);
62  double af = 0.;
63  double vf = 0.;
64  if (idAbs < 9 && idAbs%2 == 1) {
65  af = -1. + 2. * sin2tW;
66  vf = -1. + 4. * sin2tW / 3.;
67  } else if (idAbs < 9) {
68  af = 1. - 2. * sin2tW;
69  vf = 1. - 8. * sin2tW / 3.;
70  } else if (idAbs < 19 && idAbs%2 == 1) {
71  af = -1. + 2. * sin2tW;
72  vf = -1. + 4. * sin2tW;
73  }
74 
75  // Colour factor. Answer.
76  double sigma = (vf*vf + af*af) * sigma0;
77  if (idAbs < 9) sigma /= 3.;
78  return sigma;
79 
80 }
81 
82 //--------------------------------------------------------------------------
83 
84 // Select identity, colour and anticolour.
85 
86 void Sigma1ffbar2ZRight::setIdColAcol() {
87 
88  // Flavours trivial.
89  setId( id1, id2, idZR);
90 
91  // Colour flow topologies. Swap when antiquarks.
92  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
93  else setColAcol( 0, 0, 0, 0, 0, 0);
94  if (id1 < 0) swapColAcol();
95 
96 }
97 
98 //--------------------------------------------------------------------------
99 
100 // Evaluate weight for Z_R decay angle.
101 
102 double Sigma1ffbar2ZRight::weightDecay( Event& process, int iResBeg,
103  int iResEnd) {
104 
105  // Identity of mother of decaying reseonance(s).
106  int idMother = process[process[iResBeg].mother1()].idAbs();
107 
108  // For top decay hand over to standard routine.
109  if (idMother == 6)
110  return weightTopDecay( process, iResBeg, iResEnd);
111 
112  // Z_R should sit in entry 5.
113  if (iResBeg != 5 || iResEnd != 5) return 1.;
114 
115  // Couplings for in- and out-flavours.
116  double ai, vi, af, vf;
117  int idInAbs = process[3].idAbs();
118  if (idInAbs < 9 && idInAbs%2 == 1) {
119  ai = -1. + 2. * sin2tW;
120  vi = -1. + 4. * sin2tW / 3.;
121  } else if (idInAbs < 9) {
122  ai = 1. - 2. * sin2tW;
123  vi = 1. - 8. * sin2tW / 3.;
124  } else {
125  ai = -1. + 2. * sin2tW;
126  vi = -1. + 4. * sin2tW;
127  }
128  int idOutAbs = process[6].idAbs();
129  if (idOutAbs < 9 && idOutAbs%2 == 1) {
130  af = -1. + 2. * sin2tW;
131  vf = -1. + 4. * sin2tW / 3.;
132  } else if (idOutAbs < 9) {
133  af = 1. - 2. * sin2tW;
134  vf = 1. - 8. * sin2tW / 3.;
135  } else {
136  af = -1. + 2. * sin2tW;
137  vf = -1. + 4. * sin2tW;
138  }
139 
140  // Phase space factors. Reconstruct decay angle.
141  double mr1 = pow2(process[6].m()) / sH;
142  double mr2 = pow2(process[7].m()) / sH;
143  double betaf = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
144  double cosThe = (process[3].p() - process[4].p())
145  * (process[7].p() - process[6].p()) / (sH * betaf);
146 
147  // Angular weight and its maximum.
148  double wt1 = (vi*vi + ai*ai) * (vf*vf + af*af * betaf*betaf);
149  double wt2 = (1. - betaf*betaf) * (vi*vi + ai*ai) * vf*vf;
150  double wt3 = betaf * 4. * vi * ai * vf * af;
151  if (process[3].id() * process[6].id() < 0) wt3 = -wt3;
152  double wt = wt1 * (1. + cosThe*cosThe) + wt2 * (1. - cosThe*cosThe)
153  + 2. * wt3 * cosThe;
154  double wtMax = 2. * (wt1 + abs(wt3));
155 
156  // Done.
157  return wt / wtMax;
158 
159 }
160 
161 //==========================================================================
162 
163 // Sigma1ffbar2WRight class.
164 // Cross section for f fbar' -> W_R^+- (righthanded gauge boson).
165 
166 //--------------------------------------------------------------------------
167 
168 // Initialize process.
169 
170 void Sigma1ffbar2WRight::initProc() {
171 
172  // Store W_R^+- mass and width for propagator.
173  idWR = 9900024;
174  mRes = particleDataPtr->m0(idWR);
175  GammaRes = particleDataPtr->mWidth(idWR);
176  m2Res = mRes*mRes;
177  GamMRat = GammaRes / mRes;
178  thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());
179 
180  // Set pointer to particle properties and decay table.
181  particlePtr = particleDataPtr->particleDataEntryPtr(idWR);
182 
183 }
184 
185 //--------------------------------------------------------------------------
186 
187 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
188 
189 void Sigma1ffbar2WRight::sigmaKin() {
190 
191  // Common coupling factors.
192  double colQ = 3. * (1. + alpS / M_PI);
193 
194  // Reset quantities to sum. Declare variables inside loop.
195  double widOutPos = 0.;
196  double widOutNeg = 0.;
197  int id1Now, id2Now, id1Abs, id2Abs, id1Neg, id2Neg, onMode;
198  double widNow, widSecPos, widSecNeg, mf1, mf2, mr1, mr2, kinFac;
199 
200  // Loop over all W_R^+- decay channels.
201  for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
202  id1Now = particlePtr->channel(i).product(0);
203  id2Now = particlePtr->channel(i).product(1);
204  id1Abs = abs(id1Now);
205  id2Abs = abs(id2Now);
206 
207  // Check that above threshold. Phase space.
208  mf1 = particleDataPtr->m0(id1Abs);
209  mf2 = particleDataPtr->m0(id2Abs);
210  if (mH > mf1 + mf2 + MASSMARGIN) {
211  mr1 = pow2(mf1 / mH);
212  mr2 = pow2(mf2 / mH);
213  kinFac = (1. - 0.5 * (mr1 + mr2) - 0.5 * pow2(mr1 - mr2))
214  * sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 );
215 
216  // Combine kinematics with colour factor and CKM couplings.
217  widNow = kinFac;
218  if (id1Abs < 9) widNow *= colQ * couplingsPtr->V2CKMid(id1Abs, id2Abs);
219 
220  // Secondary width from top and righthanded neutrino decay.
221  id1Neg = (id1Abs < 19) ? -id1Now : id1Abs;
222  id2Neg = (id2Abs < 19) ? -id2Now : id2Abs;
223  widSecPos = particleDataPtr->resOpenFrac(id1Now, id2Now);
224  widSecNeg = particleDataPtr->resOpenFrac(id1Neg, id2Neg);
225 
226  // Add weight for channels on for all, W_R^+ and W_R^-, respectively.
227  onMode = particlePtr->channel(i).onMode();
228  if (onMode == 1 || onMode == 2) widOutPos += widNow * widSecPos;
229  if (onMode == 1 || onMode == 3) widOutNeg += widNow * widSecNeg;
230 
231  // End loop over fermions.
232  }
233  }
234 
235  // Set up Breit-Wigner. Cross section for W_R^+ and W_R^- separately.
236  double sigBW = 12. * M_PI * pow2(alpEM * thetaWRat) * sH
237  / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
238  sigma0Pos = sigBW * widOutPos;
239  sigma0Neg = sigBW * widOutNeg;
240 
241 }
242 
243 //--------------------------------------------------------------------------
244 
245 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
246 
247 double Sigma1ffbar2WRight::sigmaHat() {
248 
249  // Secondary width for W_R^+ or W_R^-. CKM and colour factors.
250  int idUp = (abs(id1)%2 == 0) ? id1 : id2;
251  double sigma = (idUp > 0) ? sigma0Pos : sigma0Neg;
252  if (abs(id1) < 9) sigma *= couplingsPtr->V2CKMid(abs(id1), abs(id2)) / 3.;
253 
254  // Answer.
255  return sigma;
256 
257 }
258 
259 //--------------------------------------------------------------------------
260 
261 // Select identity, colour and anticolour.
262 
263 void Sigma1ffbar2WRight::setIdColAcol() {
264 
265  // Sign of outgoing W_R.
266  int sign = (abs(id1)%2 == 0) ? 1 : -1;
267  if (id1 < 0) sign = -sign;
268  setId( id1, id2, idWR * sign);
269 
270  // Colour flow topologies. Swap when antiquarks.
271  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
272  else setColAcol( 0, 0, 0, 0, 0, 0);
273  if (id1 < 0) swapColAcol();
274 
275 }
276 
277 //--------------------------------------------------------------------------
278 
279 // Evaluate weight for W_R decay angle.
280 
281 double Sigma1ffbar2WRight::weightDecay( Event& process, int iResBeg,
282  int iResEnd) {
283 
284  // Identity of mother of decaying reseonance(s).
285  int idMother = process[process[iResBeg].mother1()].idAbs();
286 
287  // For top decay hand over to standard routine.
288  if (idMother == 6)
289  return weightTopDecay( process, iResBeg, iResEnd);
290 
291  // W_R should sit in entry 5.
292  if (iResBeg != 5 || iResEnd != 5) return 1.;
293 
294  // Phase space factors.
295  double mr1 = pow2(process[6].m()) / sH;
296  double mr2 = pow2(process[7].m()) / sH;
297  double betaf = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
298 
299  // Sign of asymmetry.
300  double eps = (process[3].id() * process[6].id() > 0) ? 1. : -1.;
301 
302  // Reconstruct decay angle and weight for it.
303  double cosThe = (process[3].p() - process[4].p())
304  * (process[7].p() - process[6].p()) / (sH * betaf);
305  double wtMax = 4.;
306  double wt = pow2(1. + betaf * eps * cosThe) - pow2(mr1 - mr2);
307 
308  // Done.
309  return (wt / wtMax);
310 
311 }
312 
313 //==========================================================================
314 
315 // Sigma1ll2Hchgchg class.
316 // Cross section for l l -> H_L^++-- or H_R^++-- (doubly charged Higgs).
317 
318 //--------------------------------------------------------------------------
319 
320 // Initialize process.
321 
322 void Sigma1ll2Hchgchg::initProc() {
323 
324  // Set process properties: H_L^++-- or H_R^++--.
325  if (leftRight == 1) {
326  idHLR = 9900041;
327  codeSave = 3121;
328  nameSave = "l l -> H_L^++--";
329  } else {
330  idHLR = 9900042;
331  codeSave = 3141;
332  nameSave = "l l -> H_R^++--";
333  }
334 
335  // Read in Yukawa matrix for couplings to a lepton pair.
336  yukawa[1][1] = settingsPtr->parm("LeftRightSymmmetry:coupHee");
337  yukawa[2][1] = settingsPtr->parm("LeftRightSymmmetry:coupHmue");
338  yukawa[2][2] = settingsPtr->parm("LeftRightSymmmetry:coupHmumu");
339  yukawa[3][1] = settingsPtr->parm("LeftRightSymmmetry:coupHtaue");
340  yukawa[3][2] = settingsPtr->parm("LeftRightSymmmetry:coupHtaumu");
341  yukawa[3][3] = settingsPtr->parm("LeftRightSymmmetry:coupHtautau");
342 
343  // Store H_L/R mass and width for propagator.
344  mRes = particleDataPtr->m0(idHLR);
345  GammaRes = particleDataPtr->mWidth(idHLR);
346  m2Res = mRes*mRes;
347  GamMRat = GammaRes / mRes;
348 
349  // Set pointer to particle properties and decay table.
350  particlePtr = particleDataPtr->particleDataEntryPtr(idHLR);
351 
352 }
353 
354 //--------------------------------------------------------------------------
355 
356 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
357 
358 double Sigma1ll2Hchgchg::sigmaHat() {
359 
360  // Initial state must consist of two identical-sign leptons.
361  if (id1 * id2 < 0) return 0.;
362  int id1Abs = abs(id1);
363  int id2Abs = abs(id2);
364  if (id1Abs != 11 && id1Abs != 13 && id1Abs != 15) return 0.;
365  if (id2Abs != 11 && id2Abs != 13 && id2Abs != 15) return 0.;
366 
367  // Set up Breit-Wigner, inwidth and outwidth.
368  double sigBW = 8. * M_PI / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
369  double widIn = pow2(yukawa[(id1Abs-9)/2][(id2Abs-9)/2])
370  * mH / (8. * M_PI);
371  int idSgn = (id1 < 0) ? idHLR : -idHLR;
372  double widOut = particlePtr->resWidthOpen( idSgn, mH);
373 
374  // Answer.
375  return widIn * sigBW * widOut;
376 
377 }
378 
379 //--------------------------------------------------------------------------
380 
381 // Select identity, colour and anticolour.
382 
383 void Sigma1ll2Hchgchg::setIdColAcol() {
384 
385  // Sign of outgoing H_L/R.
386  int idSgn = (id1 < 0) ? idHLR : -idHLR;
387  setId( id1, id2, idSgn);
388 
389  // No colours whatsoever.
390  setColAcol( 0, 0, 0, 0, 0, 0);
391 
392 }
393 
394 //--------------------------------------------------------------------------
395 
396 // Evaluate weight for H_L/R sequential decay angles.
397 
398 double Sigma1ll2Hchgchg::weightDecay( Event& process, int iResBeg,
399  int iResEnd) {
400 
401  // Identity of mother of decaying reseonance(s).
402  int idMother = process[process[iResBeg].mother1()].idAbs();
403 
404  // For top decay hand over to standard routine.
405  if (idMother == 6)
406  return weightTopDecay( process, iResBeg, iResEnd);
407 
408  // Else isotropic decay.
409  return 1.;
410 
411 }
412 
413 //==========================================================================
414 
415 // Sigma2lgm2Hchgchgl class.
416 // Cross section for l l -> H_L^++-- or H_R^++-- (doubly charged Higgs).
417 
418 //--------------------------------------------------------------------------
419 
420 // Initialize process.
421 
422 void Sigma2lgm2Hchgchgl::initProc() {
423 
424  // Set process properties: H_L^++-- or H_R^++-- and e/mu/tau.
425  idHLR = (leftRight == 1) ? 9900041 : 9900042;
426  codeSave = (leftRight == 1) ? 3122 : 3142;
427  if (idLep == 13) codeSave += 2;
428  if (idLep == 15) codeSave += 4;
429  if (codeSave == 3122) nameSave = "l^+- gamma -> H_L^++-- e^-+";
430  else if (codeSave == 3123) nameSave = "l^+- gamma -> H_L^++-- mu^-+";
431  else if (codeSave == 3124) nameSave = "l^+- gamma -> H_L^++-- tau^-+";
432  else if (codeSave == 3142) nameSave = "l^+- gamma -> H_R^++-- e^-+";
433  else if (codeSave == 3143) nameSave = "l^+- gamma -> H_R^++-- mu^-+";
434  else nameSave = "l^+- gamma -> H_R^++-- tau^-+";
435 
436  // Read in relevantYukawa matrix for couplings to a lepton pair.
437  if (idLep == 11) {
438  yukawa[1] = settingsPtr->parm("LeftRightSymmmetry:coupHee");
439  yukawa[2] = settingsPtr->parm("LeftRightSymmmetry:coupHmue");
440  yukawa[3] = settingsPtr->parm("LeftRightSymmmetry:coupHtaue");
441  } else if (idLep == 13) {
442  yukawa[1] = settingsPtr->parm("LeftRightSymmmetry:coupHmue");
443  yukawa[2] = settingsPtr->parm("LeftRightSymmmetry:coupHmumu");
444  yukawa[3] = settingsPtr->parm("LeftRightSymmmetry:coupHtaumu");
445  } else {
446  yukawa[1] = settingsPtr->parm("LeftRightSymmmetry:coupHtaue");
447  yukawa[2] = settingsPtr->parm("LeftRightSymmmetry:coupHtaumu");
448  yukawa[3] = settingsPtr->parm("LeftRightSymmmetry:coupHtautau");
449  }
450 
451  // Secondary open width fractions.
452  openFracPos = particleDataPtr->resOpenFrac( idHLR);
453  openFracNeg = particleDataPtr->resOpenFrac(-idHLR);
454 
455 }
456 
457 //--------------------------------------------------------------------------
458 
459 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
460 
461 double Sigma2lgm2Hchgchgl::sigmaHat() {
462 
463  // Initial state must consist of a lepton and a photon.
464  int idIn = (id2 == 22) ? id1 : id2;
465  int idInAbs = abs(idIn);
466  if (idInAbs != 11 && idInAbs != 13 && idInAbs != 15) return 0.;
467 
468  // Incoming squared lepton mass.
469  double s1 = pow2( particleDataPtr->m0(idInAbs) );
470 
471  // Kinematical expressions.
472  double smm1 = 8. * (sH + tH - s3) * (sH + tH - 2. * s3 - s1 - s4)
473  / pow2(uH - s3);
474  double smm2 = 2. * ( (2. * s3 - 3. * s1) * s4 + (s1 - 2. * s4) * tH
475  - (tH - s4) * sH ) / pow2(tH - s4);
476  double smm3 = 2. * ( (2. * s3 - 3. * s4 + tH) * s1
477  - (2. * s1 - s4 + tH) * sH ) / pow2(sH - s1);
478  double smm12 = 4. * ( (2. * s1 - s4 - 2. * s3 + tH) * sH
479  + (tH - 3. * s3 - 3. * s4) * tH + (2. * s3 - 2. * s1
480  + 3. * s4) * s3 ) / ( (uH - s3) * (tH - s4) );
481  double smm13 = -4. * ( (tH + s1 - 2. * s4) * tH - (s3 + 3. * s1 - 2. * s4)
482  * s3 + (s3 + 3. * s1 + tH) * sH - pow2(tH - s3 + sH) )
483  / ( (uH - s3) * (sH - s1) );
484  double smm23 = -4. * ( (s1 - s4 + s3) * tH - s3*s3 + s3 * (s1 + s4)
485  - 3. * s1 * s4 - (s1 - s4 - s3 + tH) * sH)
486  / ( (sH - s1) * (tH - s4) );
487  double sigma = alpEM * pow2(sH / (sH - s1) ) * (smm1 + smm2 + smm3
488  + smm12 + smm13 + smm23) / (4. * sH2);
489 
490  // Lepton Yukawa and secondary widths.
491  sigma *= pow2(yukawa[(idInAbs-9)/2]);
492  sigma *= (idIn < 0) ? openFracPos : openFracNeg;
493 
494  // Answer.
495  return sigma;
496 
497 }
498 
499 //--------------------------------------------------------------------------
500 
501 // Select identity, colour and anticolour.
502 
503 void Sigma2lgm2Hchgchgl::setIdColAcol() {
504 
505  // Sign of outgoing H_L/R.
506  int idIn = (id2 == 22) ? id1 : id2;
507  int idSgn = (idIn < 0) ? idHLR : -idHLR;
508  int idOut = (idIn < 0) ? idLep : -idLep;
509  setId( id1, id2, idSgn, idOut);
510 
511  // tHat is defined between incoming lepton and outgoing Higgs.
512  if (id1 == 22) swapTU = true;
513 
514  // No colours whatsoever.
515  setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
516 
517 }
518 
519 //--------------------------------------------------------------------------
520 
521 // Evaluate weight for H_L/R sequential decay angles.
522 
523 double Sigma2lgm2Hchgchgl::weightDecay( Event& process, int iResBeg,
524  int iResEnd) {
525 
526  // Identity of mother of decaying reseonance(s).
527  int idMother = process[process[iResBeg].mother1()].idAbs();
528 
529  // For top decay hand over to standard routine.
530  if (idMother == 6)
531  return weightTopDecay( process, iResBeg, iResEnd);
532 
533  // Else isotropic decay.
534  return 1.;
535 
536 }
537 
538 //==========================================================================
539 
540 // Sigma3ff2HchgchgfftWW class.
541 // Cross section for f_1 f_2 -> H_(L/R)^++-- f_3 f_4 (W+- W+- fusion).
542 
543 //--------------------------------------------------------------------------
544 
545 // Initialize process.
546 
547 void Sigma3ff2HchgchgfftWW::initProc() {
548 
549  // Set process properties: H_L^++-- or H_R^++--.
550  if (leftRight == 1) {
551  idHLR = 9900041;
552  codeSave = 3125;
553  nameSave = "f_1 f_2 -> H_L^++-- f_3 f_4 (W+- W+- fusion)";
554  } else {
555  idHLR = 9900042;
556  codeSave = 3145;
557  nameSave = "f_1 f_2 -> H_R^++-- f_3 f_4 (W+- W+- fusion)";
558  }
559 
560  // Common fixed mass and coupling factor.
561  double mW = particleDataPtr->m0(24);
562  double mWR = particleDataPtr->m0(9900024);
563  mWS = (leftRight == 1) ? pow2(mW) : pow2(mWR);
564  double gL = settingsPtr->parm("LeftRightSymmmetry:gL");
565  double gR = settingsPtr->parm("LeftRightSymmmetry:gR");
566  double vL = settingsPtr->parm("LeftRightSymmmetry:vL");
567  prefac = (leftRight == 1) ? pow2(pow4(gL) * vL)
568  : 2. * pow2(pow3(gR) * mWR);
569  // Secondary open width fractions.
570  openFracPos = particleDataPtr->resOpenFrac( idHLR);
571  openFracNeg = particleDataPtr->resOpenFrac(-idHLR);
572 
573 }
574 
575 //--------------------------------------------------------------------------
576 
577 // Evaluate sigmaHat(sHat), part independent of incoming flavour.
578 
579 void Sigma3ff2HchgchgfftWW::sigmaKin() {
580 
581  // Required four-vector products.
582  double pp12 = 0.5 * sH;
583  double pp14 = 0.5 * mH * p4cm.pNeg();
584  double pp15 = 0.5 * mH * p5cm.pNeg();
585  double pp24 = 0.5 * mH * p4cm.pPos();
586  double pp25 = 0.5 * mH * p5cm.pPos();
587  double pp45 = p4cm * p5cm;
588 
589  // Cross section: kinematics part. Combine with couplings.
590  double propT = 1. / ( (2. * pp14 + mWS) * (2. * pp25 + mWS) );
591  double propU = 1. / ( (2. * pp24 + mWS) * (2. * pp15 + mWS) );
592  sigma0TU = prefac * pp12 * pp45 * pow2(propT + propU);
593  sigma0T = prefac * pp12 * pp45 * 2. * pow2(propT);
594 
595 }
596 
597 //--------------------------------------------------------------------------
598 
599 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
600 
601 double Sigma3ff2HchgchgfftWW::sigmaHat() {
602 
603  // Do not allow creation of righthanded neutrinos for H_R.
604  int id1Abs = abs(id1);
605  int id2Abs = abs(id2);
606  if ( leftRight == 2 && (id1Abs > 10 || id2Abs > 10) ) return 0.;
607 
608  // Many flavour combinations not possible because of charge.
609  int chg1 = (( id1Abs%2 == 0 && id1 > 0)
610  || (id1Abs%2 == 1 && id1 < 0) ) ? 1 : -1;
611  int chg2 = (( id2Abs%2 == 0 && id2 > 0)
612  || (id2Abs%2 == 1 && id2 < 0) ) ? 1 : -1;
613  if (abs(chg1 + chg2) != 2) return 0.;
614 
615  // Basic cross section. CKM factors for final states.
616  double sigma = (id2 == id1 && id1Abs > 10) ? sigma0TU : sigma0T;
617  sigma *= couplingsPtr->V2CKMsum(id1Abs)
618  * couplingsPtr->V2CKMsum(id2Abs);
619 
620  // Secondary width for H0.
621  sigma *= (chg1 + chg2 == 2) ? openFracPos : openFracNeg;
622 
623  // Spin-state extra factor 2 per incoming neutrino.
624  if (id1Abs == 12 || id1Abs == 14 || id1Abs == 16) sigma *= 2.;
625  if (id2Abs == 12 || id2Abs == 14 || id2Abs == 16) sigma *= 2.;
626 
627  // Answer.
628  return sigma;
629 
630 }
631 
632 //--------------------------------------------------------------------------
633 
634 // Select identity, colour and anticolour.
635 
636 void Sigma3ff2HchgchgfftWW::setIdColAcol() {
637 
638  // Pick out-flavours by relative CKM weights.
639  int id1Abs = abs(id1);
640  int id2Abs = abs(id2);
641  id4 = couplingsPtr->V2CKMpick(id1);
642  id5 = couplingsPtr->V2CKMpick(id2);
643 
644  // Find charge of Higgs .
645  id3 = (( id1Abs%2 == 0 && id1 > 0) || (id1Abs%2 == 1 && id1 < 0) )
646  ? idHLR : -idHLR;
647  setId( id1, id2, id3, id4, id5);
648 
649  // Colour flow topologies. Swap when antiquarks.
650  if (id1Abs < 9 && id2Abs < 9 && id1*id2 > 0)
651  setColAcol( 1, 0, 2, 0, 0, 0, 1, 0, 2, 0);
652  else if (id1Abs < 9 && id2Abs < 9)
653  setColAcol( 1, 0, 0, 2, 0, 0, 1, 0, 0, 2);
654  else if (id1Abs < 9) setColAcol( 1, 0, 0, 0, 0, 0, 1, 0, 0, 0);
655  else if (id2Abs < 9) setColAcol( 0, 0, 1, 0, 0, 0, 0, 0, 1, 0);
656  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
657  if ( (id1Abs < 9 && id1 < 0) || (id1Abs > 10 && id2 < 0) )
658  swapColAcol();
659 
660 }
661 
662 //--------------------------------------------------------------------------
663 
664 // Evaluate weight for decay angles.
665 
666 double Sigma3ff2HchgchgfftWW::weightDecay( Event& process, int iResBeg,
667  int iResEnd) {
668 
669  // Identity of mother of decaying reseonance(s).
670  int idMother = process[process[iResBeg].mother1()].idAbs();
671 
672  // For top decay hand over to standard routine.
673  if (idMother == 6)
674  return weightTopDecay( process, iResBeg, iResEnd);
675 
676  // Else done.
677  return 1.;
678 
679 }
680 
681 //==========================================================================
682 
683 // Sigma2ffbar2HchgchgHchgchg class.
684 // Cross section for f fbar -> H_(L/R)^++ H_(L/R)^-- (doubly charged Higgs).
685 
686 //--------------------------------------------------------------------------
687 
688 // Initialize process.
689 
690 void Sigma2ffbar2HchgchgHchgchg::initProc() {
691 
692  // Set process properties: H_L^++ H_L^-- or H_R^++ H_R^--.
693  if (leftRight == 1) {
694  idHLR = 9900041;
695  codeSave = 3126;
696  nameSave = "f fbar -> H_L^++ H_L^--";
697  } else {
698  idHLR = 9900042;
699  codeSave = 3146;
700  nameSave = "f fbar -> H_R^++ H_R^--";
701  }
702 
703  // Read in Yukawa matrix for couplings to a lepton pair.
704  yukawa[1][1] = settingsPtr->parm("LeftRightSymmmetry:coupHee");
705  yukawa[2][1] = settingsPtr->parm("LeftRightSymmmetry:coupHmue");
706  yukawa[2][2] = settingsPtr->parm("LeftRightSymmmetry:coupHmumu");
707  yukawa[3][1] = settingsPtr->parm("LeftRightSymmmetry:coupHtaue");
708  yukawa[3][2] = settingsPtr->parm("LeftRightSymmmetry:coupHtaumu");
709  yukawa[3][3] = settingsPtr->parm("LeftRightSymmmetry:coupHtautau");
710 
711  // Electroweak parameters.
712  mRes = particleDataPtr->m0(23);
713  GammaRes = particleDataPtr->mWidth(23);
714  m2Res = mRes*mRes;
715  GamMRat = GammaRes / mRes;
716  sin2tW = couplingsPtr->sin2thetaW();
717  preFac = (1. - 2. * sin2tW) / ( 8. * sin2tW * (1. - sin2tW) );
718 
719  // Open fraction from secondary widths.
720  openFrac = particleDataPtr->resOpenFrac( idHLR, -idHLR);
721 
722 }
723 
724 //--------------------------------------------------------------------------
725 
726 // Evaluate sigmaHat(sHat), including incoming flavour dependence.
727 
728 double Sigma2ffbar2HchgchgHchgchg::sigmaHat() {
729 
730  // Electroweak couplings to gamma^*/Z^0.
731  int idAbs = abs(id1);
732  double ei = couplingsPtr->ef(idAbs);
733  double vi = couplingsPtr->vf(idAbs);
734  double ai = couplingsPtr->af(idAbs);
735 
736  // Part via gamma^*/Z^0 propagator.No Z^0 coupling to H_R.
737  double resProp = 1. / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
738  double sigma = 8. * pow2(alpEM) * ei*ei / sH2;
739  if (leftRight == 1) sigma += 8. * pow2(alpEM)
740  * (2. * ei * vi * preFac * (sH - m2Res) * resProp / sH
741  + (vi * vi + ai * ai) * pow2(preFac) * resProp);
742 
743  // Part via t-channel lepton + interference; sum over possibilities.
744  if (idAbs == 11 || idAbs == 13 || idAbs == 15) {
745  double yuk2Sum;
746  if (idAbs == 11) yuk2Sum
747  = pow2(yukawa[1][1]) + pow2(yukawa[2][1]) + pow2(yukawa[3][1]);
748  else if (idAbs == 13) yuk2Sum
749  = pow2(yukawa[2][1]) + pow2(yukawa[2][2]) + pow2(yukawa[3][2]);
750  else yuk2Sum
751  = pow2(yukawa[3][1]) + pow2(yukawa[3][2]) + pow2(yukawa[3][3]);
752  yuk2Sum /= 4. * M_PI;
753  sigma += 8. * alpEM * ei * yuk2Sum / (sH * tH)
754  + 4. * pow2(yuk2Sum) / tH2;
755  if (leftRight == 1) sigma += 8. * alpEM * (vi + ai) * yuk2Sum
756  * preFac * (sH - m2Res) * resProp / tH;
757  }
758 
759  // Common kinematical factor. Colour factor.
760  sigma *= M_PI * (tH * uH - s3 * s4) / sH2;
761  if (idAbs < 9) sigma /= 3.;
762 
763  // Answer.
764  return sigma;
765 
766 }
767 
768 //--------------------------------------------------------------------------
769 
770 // Select identity, colour and anticolour.
771 
772 void Sigma2ffbar2HchgchgHchgchg::setIdColAcol() {
773 
774  // Outgoing flavours trivial.
775  setId( id1, id2, idHLR, -idHLR);
776 
777  // tHat is defined between incoming fermion and outgoing H--.
778  if (id1 > 0) swapTU = true;
779 
780  // No colours at all or one flow topology. Swap if first is antiquark.
781  if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0, 0, 0);
782  else setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
783  if (id1 < 0) swapColAcol();
784 
785 }
786 
787 //--------------------------------------------------------------------------
788 
789 // Evaluate weight for H_L/R sequential decay angles.
790 
791 double Sigma2ffbar2HchgchgHchgchg::weightDecay( Event& process,
792  int iResBeg, int iResEnd) {
793 
794  // Identity of mother of decaying reseonance(s).
795  int idMother = process[process[iResBeg].mother1()].idAbs();
796 
797  // For top decay hand over to standard routine.
798  if (idMother == 6)
799  return weightTopDecay( process, iResBeg, iResEnd);
800 
801  // Else isotropic decay.
802  return 1.;
803 
804 }
805 
806 //==========================================================================
807 
808 } // end namespace Pythia8
Event
Definition: AgUStep.h:26
Generated by   doxygen 1.8.5