← Back to branch summary

~ubuntu-branches/ubuntu/trusty/pythia8/trusty-proposed

« back to all changes in this revision

Viewing changes to src/SigmaLeftRightSym.cc

  • Committer: Package Import Robot
  • Author(s): Lifeng Sun
  • Date: 2012-05-22 11:43:00 UTC
  • Revision ID: package-import@ubuntu.com-20120522114300-0jvsv2vl4o2bo435
Tags: upstream-8.1.65
Import upstream version 8.1.65

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/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