using option '-treename HLTesdTree' for EsdCollector, adding default parameter for...
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8130 / src / SigmaLeftRightSym.cxx
CommitLineData
5ad4eb21 1// SigmaLeftRightSym.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// left-right-symmetry simulation classes.
8
9#include "SigmaLeftRightSym.h"
10
11namespace 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
22void Sigma1ffbar2ZRight::initProc() {
23
24 // Store Z_R mass and width for propagator.
25 idZR = 9900023;
26 mRes = ParticleDataTable::m0(idZR);
27 GammaRes = ParticleDataTable::mWidth(idZR);
28 m2Res = mRes*mRes;
29 GamMRat = GammaRes / mRes;
30 sin2tW = CoupEW::sin2thetaW();
31
32 // Set pointer to particle properties and decay table.
33 ZRPtr = ParticleDataTable::particleDataPtr(idZR);
34
35}
36
37//*********
38
39// Evaluate sigmaHat(sHat), part independent of incoming flavour.
40
41void 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
58double 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
86void 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
102double 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
170void Sigma1ffbar2WRight::initProc() {
171
172 // Store W_R^+- mass and width for propagator.
173 idWR = 9900024;
174 mRes = ParticleDataTable::m0(idWR);
175 GammaRes = ParticleDataTable::mWidth(idWR);
176 m2Res = mRes*mRes;
177 GamMRat = GammaRes / mRes;
178 thetaWRat = 1. / (12. * CoupEW::sin2thetaW());
179
180 // Set pointer to particle properties and decay table.
181 particlePtr = ParticleDataTable::particleDataPtr(idWR);
182
183}
184
185//*********
186
187// Evaluate sigmaHat(sHat), part independent of incoming flavour.
188
189void 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->decay.size(); ++i) {
202 widNow = 0.;
203 id1Now = particlePtr->decay[i].product(0);
204 id2Now = particlePtr->decay[i].product(1);
205 id1Abs = abs(id1Now);
206 id2Abs = abs(id2Now);
207
208 // Check that above threshold. Phase space.
209 mf1 = ParticleDataTable::m0(id1Abs);
210 mf2 = ParticleDataTable::m0(id2Abs);
211 if (mH > mf1 + mf2 + MASSMARGIN) {
212 mr1 = pow2(mf1 / mH);
213 mr2 = pow2(mf2 / mH);
214 kinFac = (1. - 0.5 * (mr1 + mr2) - 0.5 * pow2(mr1 - mr2))
215 * sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 );
216
217 // Combine kinematics with colour factor and CKM couplings.
218 widNow = kinFac;
219 if (id1Abs < 9) widNow *= colQ * VCKM::V2id(id1Abs, id2Abs);
220
221 // Secondary width from top and righthanded neutrino decay.
222 id1Neg = (id1Abs < 19) ? -id1Now : id1Abs;
223 id2Neg = (id2Abs < 19) ? -id2Now : id2Abs;
224 widSecPos = ParticleDataTable::resOpenFrac(id1Now, id2Now);
225 widSecNeg = ParticleDataTable::resOpenFrac(id1Neg, id2Neg);
226
227 // Add weight for channels on for all, W_R^+ and W_R^-, respectively.
228 onMode = particlePtr->decay[i].onMode();
229 if (onMode == 1 || onMode == 2) widOutPos += widNow * widSecPos;
230 if (onMode == 1 || onMode == 3) widOutNeg += widNow * widSecNeg;
231
232 // End loop over fermions.
233 }
234 }
235
236 // Set up Breit-Wigner. Cross section for W_R^+ and W_R^- separately.
237 double sigBW = 12. * M_PI * pow2(alpEM * thetaWRat) * sH
238 / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
239 sigma0Pos = sigBW * widOutPos;
240 sigma0Neg = sigBW * widOutNeg;
241
242}
243
244//*********
245
246// Evaluate sigmaHat(sHat), including incoming flavour dependence.
247
248double Sigma1ffbar2WRight::sigmaHat() {
249
250 // Secondary width for W_R^+ or W_R^-. CKM and colour factors.
251 int idUp = (abs(id1)%2 == 0) ? id1 : id2;
252 double sigma = (idUp > 0) ? sigma0Pos : sigma0Neg;
253 if (abs(id1) < 9) sigma *= VCKM::V2id(abs(id1), abs(id2)) / 3.;
254
255 // Answer.
256 return sigma;
257
258}
259
260//*********
261
262// Select identity, colour and anticolour.
263
264void Sigma1ffbar2WRight::setIdColAcol() {
265
266 // Sign of outgoing W_R.
267 int sign = (abs(id1)%2 == 0) ? 1 : -1;
268 if (id1 < 0) sign = -sign;
269 setId( id1, id2, idWR * sign);
270
271 // Colour flow topologies. Swap when antiquarks.
272 if (abs(id1) < 9) setColAcol( 1, 0, 0, 1, 0, 0);
273 else setColAcol( 0, 0, 0, 0, 0, 0);
274 if (id1 < 0) swapColAcol();
275
276}
277
278//*********
279
280// Evaluate weight for W_R decay angle.
281
282double Sigma1ffbar2WRight::weightDecay( Event& process, int iResBeg,
283 int iResEnd) {
284
285 // Identity of mother of decaying reseonance(s).
286 int idMother = process[process[iResBeg].mother1()].idAbs();
287
288 // For top decay hand over to standard routine.
289 if (idMother == 6)
290 return weightTopDecay( process, iResBeg, iResEnd);
291
292 // W_R should sit in entry 5.
293 if (iResBeg != 5 || iResEnd != 5) return 1.;
294
295 // Phase space factors.
296 double mr1 = pow2(process[6].m()) / sH;
297 double mr2 = pow2(process[7].m()) / sH;
298 double betaf = sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2);
299
300 // Sign of asymmetry.
301 double eps = (process[3].id() * process[6].id() > 0) ? 1. : -1.;
302
303 // Reconstruct decay angle and weight for it.
304 double cosThe = (process[3].p() - process[4].p())
305 * (process[7].p() - process[6].p()) / (sH * betaf);
306 double wtMax = 4.;
307 double wt = pow2(1. + betaf * eps * cosThe) - pow2(mr1 - mr2);
308
309 // Done.
310 return (wt / wtMax);
311
312}
313
314//**************************************************************************
315
316// Sigma1ll2Hchgchg class.
317// Cross section for l l -> H_L^++-- or H_R^++-- (doubly charged Higgs).
318
319//*********
320
321// Initialize process.
322
323void Sigma1ll2Hchgchg::initProc() {
324
325 // Set process properties: H_L^++-- or H_R^++--.
326 if (leftRight == 1) {
327 idHLR = 9900041;
328 codeSave = 3121;
329 nameSave = "l l -> H_L^++--";
330 } else {
331 idHLR = 9900042;
332 codeSave = 3141;
333 nameSave = "l l -> H_R^++--";
334 }
335
336 // Read in Yukawa matrix for couplings to a lepton pair.
337 yukawa[1][1] = Settings::parm("LeftRightSymmmetry:coupHee");
338 yukawa[2][1] = Settings::parm("LeftRightSymmmetry:coupHmue");
339 yukawa[2][2] = Settings::parm("LeftRightSymmmetry:coupHmumu");
340 yukawa[3][1] = Settings::parm("LeftRightSymmmetry:coupHtaue");
341 yukawa[3][2] = Settings::parm("LeftRightSymmmetry:coupHtaumu");
342 yukawa[3][3] = Settings::parm("LeftRightSymmmetry:coupHtautau");
343
344 // Store H_L/R mass and width for propagator.
345 mRes = ParticleDataTable::m0(idHLR);
346 GammaRes = ParticleDataTable::mWidth(idHLR);
347 m2Res = mRes*mRes;
348 GamMRat = GammaRes / mRes;
349
350 // Set pointer to particle properties and decay table.
351 particlePtr = ParticleDataTable::particleDataPtr(idHLR);
352
353}
354
355//*********
356
357// Evaluate sigmaHat(sHat), including incoming flavour dependence.
358
359double Sigma1ll2Hchgchg::sigmaHat() {
360
361 // Initial state must consist of two identical-sign leptons.
362 if (id1 * id2 < 0) return 0.;
363 int id1Abs = abs(id1);
364 int id2Abs = abs(id2);
365 if (id1Abs != 11 && id1Abs != 13 && id1Abs != 15) return 0.;
366 if (id2Abs != 11 && id2Abs != 13 && id2Abs != 15) return 0.;
367
368 // Set up Breit-Wigner, inwidth and outwidth.
369 double sigBW = 8. * M_PI / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
370 double widIn = pow2(yukawa[(id1Abs-9)/2][(id2Abs-9)/2])
371 * mH / (8. * M_PI);
372 int idSgn = (id1 < 0) ? idHLR : -idHLR;
373 double widOut = particlePtr->resWidthOpen( idSgn, mH);
374
375 // Answer.
376 return widIn * sigBW * widOut;
377
378}
379
380//*********
381
382// Select identity, colour and anticolour.
383
384void Sigma1ll2Hchgchg::setIdColAcol() {
385
386 // Sign of outgoing H_L/R.
387 int idSgn = (id1 < 0) ? idHLR : -idHLR;
388 setId( id1, id2, idSgn);
389
390 // No colours whatsoever.
391 setColAcol( 0, 0, 0, 0, 0, 0);
392
393}
394
395//*********
396
397// Evaluate weight for H_L/R sequential decay angles.
398
399double Sigma1ll2Hchgchg::weightDecay( Event& process, int iResBeg,
400 int iResEnd) {
401
402 // Identity of mother of decaying reseonance(s).
403 int idMother = process[process[iResBeg].mother1()].idAbs();
404
405 // For top decay hand over to standard routine.
406 if (idMother == 6)
407 return weightTopDecay( process, iResBeg, iResEnd);
408
409 // Else isotropic decay.
410 return 1.;
411
412}
413
414//**************************************************************************
415
416// Sigma2lgm2Hchgchgl class.
417// Cross section for l l -> H_L^++-- or H_R^++-- (doubly charged Higgs).
418
419//*********
420
421// Initialize process.
422
423void Sigma2lgm2Hchgchgl::initProc() {
424
425 // Set process properties: H_L^++-- or H_R^++-- and e/mu/tau.
426 idHLR = (leftRight == 1) ? 9900041 : 9900042;
427 codeSave = (leftRight == 1) ? 3122 : 3142;
428 if (idLep == 13) codeSave += 2;
429 if (idLep == 15) codeSave += 4;
430 if (codeSave == 3122) nameSave = "l^+- gamma -> H_L^++-- e^-+";
431 else if (codeSave == 3123) nameSave = "l^+- gamma -> H_L^++-- mu^-+";
432 else if (codeSave == 3124) nameSave = "l^+- gamma -> H_L^++-- tau^-+";
433 else if (codeSave == 3142) nameSave = "l^+- gamma -> H_R^++-- e^-+";
434 else if (codeSave == 3143) nameSave = "l^+- gamma -> H_R^++-- mu^-+";
435 else nameSave = "l^+- gamma -> H_R^++-- tau^-+";
436
437 // Read in relevantYukawa matrix for couplings to a lepton pair.
438 if (idLep == 11) {
439 yukawa[1] = Settings::parm("LeftRightSymmmetry:coupHee");
440 yukawa[2] = Settings::parm("LeftRightSymmmetry:coupHmue");
441 yukawa[3] = Settings::parm("LeftRightSymmmetry:coupHtaue");
442 } else if (idLep == 13) {
443 yukawa[1] = Settings::parm("LeftRightSymmmetry:coupHmue");
444 yukawa[2] = Settings::parm("LeftRightSymmmetry:coupHmumu");
445 yukawa[3] = Settings::parm("LeftRightSymmmetry:coupHtaumu");
446 } else {
447 yukawa[1] = Settings::parm("LeftRightSymmmetry:coupHtaue");
448 yukawa[2] = Settings::parm("LeftRightSymmmetry:coupHtaumu");
449 yukawa[3] = Settings::parm("LeftRightSymmmetry:coupHtautau");
450 }
451
452 // Secondary open width fractions.
453 openFracPos = ParticleDataTable::resOpenFrac( idHLR);
454 openFracNeg = ParticleDataTable::resOpenFrac(-idHLR);
455
456}
457
458//*********
459
460// Evaluate sigmaHat(sHat), including incoming flavour dependence.
461
462double Sigma2lgm2Hchgchgl::sigmaHat() {
463
464 // Initial state must consist of a lepton and a photon.
465 int idIn = (id2 == 22) ? id1 : id2;
466 int idInAbs = abs(idIn);
467 if (idInAbs != 11 && idInAbs != 13 && idInAbs != 15) return 0.;
468
469 // Incoming squared lepton mass.
470 double s1 = pow2( ParticleDataTable::m0(idInAbs) );
471
472 // Kinematical expressions.
473 double smm1 = 8. * (sH + tH - s3) * (sH + tH - 2. * s3 - s1 - s4)
474 / pow2(uH - s3);
475 double smm2 = 2. * ( (2. * s3 - 3. * s1) * s4 + (s1 - 2. * s4) * tH
476 - (tH - s4) * sH ) / pow2(tH - s4);
477 double smm3 = 2. * ( (2. * s3 - 3. * s4 + tH) * s1
478 - (2. * s1 - s4 + tH) * sH ) / pow2(sH - s1);
479 double smm12 = 4. * ( (2. * s1 - s4 - 2. * s3 + tH) * sH
480 + (tH - 3. * s3 - 3. * s4) * tH + (2. * s3 - 2. * s1
481 + 3. * s4) * s3 ) / ( (uH - s3) * (tH - s4) );
482 double smm13 = -4. * ( (tH + s1 - 2. * s4) * tH - (s3 + 3. * s1 - 2. * s4)
483 * s3 + (s3 + 3. * s1 + tH) * sH - pow2(tH - s3 + sH) )
484 / ( (uH - s3) * (sH - s1) );
485 double smm23 = -4. * ( (s1 - s4 + s3) * tH - s3*s3 + s3 * (s1 + s4)
486 - 3. * s1 * s4 - (s1 - s4 - s3 + tH) * sH)
487 / ( (sH - s1) * (tH - s4) );
488 double sigma = alpEM * pow2(sH / (sH - s1) ) * (smm1 + smm2 + smm3
489 + smm12 + smm13 + smm23) / (4. * sH2);
490
491 // Lepton Yukawa and secondary widths.
492 sigma *= pow2(yukawa[(idInAbs-9)/2]);
493 sigma *= (idIn < 0) ? openFracPos : openFracNeg;
494
495 // Answer.
496 return sigma;
497
498}
499
500//*********
501
502// Select identity, colour and anticolour.
503
504void Sigma2lgm2Hchgchgl::setIdColAcol() {
505
506 // Sign of outgoing H_L/R.
507 int idIn = (id2 == 22) ? id1 : id2;
508 int idSgn = (idIn < 0) ? idHLR : -idHLR;
509 int idOut = (idIn < 0) ? idLep : -idLep;
510 setId( id1, id2, idSgn, idOut);
511
512 // tHat is defined between incoming lepton and outgoing Higgs.
513 if (id1 == 22) swapTU = true;
514
515 // No colours whatsoever.
516 setColAcol( 0, 0, 0, 0, 0, 0, 0, 0);
517
518}
519
520//*********
521
522// Evaluate weight for H_L/R sequential decay angles.
523
524double Sigma2lgm2Hchgchgl::weightDecay( Event& process, int iResBeg,
525 int iResEnd) {
526
527 // Identity of mother of decaying reseonance(s).
528 int idMother = process[process[iResBeg].mother1()].idAbs();
529
530 // For top decay hand over to standard routine.
531 if (idMother == 6)
532 return weightTopDecay( process, iResBeg, iResEnd);
533
534 // Else isotropic decay.
535 return 1.;
536
537}
538
539//**************************************************************************
540
541// Sigma3ff2HchgchgfftWW class.
542// Cross section for f_1 f_2 -> H0 f_3 f_4 (W+ W- fusion of SM Higgs).
543
544//*********
545
546// Initialize process.
547
548void Sigma3ff2HchgchgfftWW::initProc() {
549
550 // Set process properties: H_L^++-- or H_R^++--.
551 if (leftRight == 1) {
552 idHLR = 9900041;
553 codeSave = 3125;
554 nameSave = "f_1 f_2 -> H_L^++-- f_3 f_4 (W+- W+- fusion)";
555 } else {
556 idHLR = 9900042;
557 codeSave = 3145;
558 nameSave = "f_1 f_2 -> H_R^++-- f_3 f_4 (W+- W+- fusion)";
559 }
560
561 // Common fixed mass and coupling factor.
562 double mW = ParticleDataTable::m0(24);
563 double mWR = ParticleDataTable::m0(9900024);
564 mWS = (leftRight == 1) ? pow2(mW) : pow2(mWR);
565 double gL = Settings::parm("LeftRightSymmmetry:gL");
566 double gR = Settings::parm("LeftRightSymmmetry:gR");
567 double vL = Settings::parm("LeftRightSymmmetry:vL");
568 prefac = (leftRight == 1) ? pow2(pow4(gL) * vL)
569 : 2. * pow2(pow3(gR) * mWR);
570 // Secondary open width fractions.
571 openFracPos = ParticleDataTable::resOpenFrac( idHLR);
572 openFracNeg = ParticleDataTable::resOpenFrac(-idHLR);
573
574}
575
576//*********
577
578// Evaluate sigmaHat(sHat), part independent of incoming flavour.
579
580void Sigma3ff2HchgchgfftWW::sigmaKin() {
581
582 // Required four-vector products.
583 double pp12 = 0.5 * sH;
584 double pp14 = 0.5 * mH * p4cm.pMinus();
585 double pp15 = 0.5 * mH * p5cm.pMinus();
586 double pp24 = 0.5 * mH * p5cm.pPlus();
587 double pp25 = 0.5 * mH * p5cm.pPlus();
588 double pp45 = p4cm * p5cm;
589
590 // Cross section: kinematics part. Combine with couplings.
591 double propT = 1. / ( (2. * pp14 + mWS) * (2. * pp25 + mWS) );
592 double propU = 1. / ( (2. * pp24 + mWS) * (2. * pp15 + mWS) );
593 sigma0TU = prefac * pp12 * pp45 * pow2(propT + propU);
594 sigma0T = prefac * pp12 * pp45 * 2. * pow2(propT);
595
596}
597
598//*********
599
600// Evaluate sigmaHat(sHat), including incoming flavour dependence.
601
602double Sigma3ff2HchgchgfftWW::sigmaHat() {
603
604 // Do not allow creation of righthanded neutrinos for H_R.
605 int id1Abs = abs(id1);
606 int id2Abs = abs(id2);
607 if ( leftRight == 2 && (id1Abs > 10 || id2Abs > 10) ) return 0.;
608
609 // Many flavour combinations not possible because of charge.
610 int chg1 = (( id1Abs%2 == 0 && id1 > 0)
611 || (id1Abs%2 == 1 && id1 < 0) ) ? 1 : -1;
612 int chg2 = (( id2Abs%2 == 0 && id2 > 0)
613 || (id2Abs%2 == 1 && id2 < 0) ) ? 1 : -1;
614 if (abs(chg1 + chg2) != 2) return 0.;
615
616 // Basic cross section. CKM factors for final states.
617 double sigma = (id2 == id1 && id1Abs > 10) ? sigma0TU : sigma0T;
618 sigma *= VCKM::V2sum(id1Abs) * VCKM::V2sum(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
636void Sigma3ff2HchgchgfftWW::setIdColAcol() {
637
638 // Pick out-flavours by relative CKM weights.
639 int id1Abs = abs(id1);
640 int id2Abs = abs(id2);
641 id4 = VCKM::V2pick(id1);
642 id5 = VCKM::V2pick(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
666double 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
690void 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] = Settings::parm("LeftRightSymmmetry:coupHee");
705 yukawa[2][1] = Settings::parm("LeftRightSymmmetry:coupHmue");
706 yukawa[2][2] = Settings::parm("LeftRightSymmmetry:coupHmumu");
707 yukawa[3][1] = Settings::parm("LeftRightSymmmetry:coupHtaue");
708 yukawa[3][2] = Settings::parm("LeftRightSymmmetry:coupHtaumu");
709 yukawa[3][3] = Settings::parm("LeftRightSymmmetry:coupHtautau");
710
711 // Electroweak parameters.
712 mRes = ParticleDataTable::m0(23);
713 GammaRes = ParticleDataTable::mWidth(23);
714 m2Res = mRes*mRes;
715 GamMRat = GammaRes / mRes;
716 sin2tW = CoupEW::sin2thetaW();
717 preFac = (1. - 2. * sin2tW) / ( 8. * sin2tW * (1. - sin2tW) );
718
719 // Open fraction from secondary widths.
720 openFrac = ParticleDataTable::resOpenFrac( idHLR, -idHLR);
721
722}
723
724//*********
725
726// Evaluate sigmaHat(sHat), including incoming flavour dependence.
727
728double Sigma2ffbar2HchgchgHchgchg::sigmaHat() {
729
730 // Electroweak couplings to gamma^*/Z^0.
731 int idAbs = abs(id1);
732 double ei = CoupEW::ef(idAbs);
733 double vi = CoupEW::vf(idAbs);
734 double ai = CoupEW::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
772void 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
791double 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