StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
ResonanceWidths.cc
1 // ResonanceWidths.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
7 // the ResonanceWidths class and classes derived from it.
8 
9 #include "ResonanceWidths.h"
10 #include "PythiaComplex.h"
11 
12 namespace Pythia8 {
13 
14 //==========================================================================
15 
16 // The ResonanceWidths class.
17 // Base class for the various resonances.
18 
19 //--------------------------------------------------------------------------
20 
21 // Constants: could be changed here if desired, but normally should not.
22 // These are of technical nature, as described for each.
23 
24 // Number of points in integration direction for numInt routines.
25 const int ResonanceWidths::NPOINT = 100;
26 
27 // The sum of product masses must not be too close to the resonance mass.
28 const double ResonanceWidths::MASSMARGIN = 0.1;
29 
30 //--------------------------------------------------------------------------
31 
32 // Initialize data members.
33 // Calculate and store partial and total widths at the nominal mass.
34 
35 bool ResonanceWidths::init(Info* infoPtrIn, Settings* settingsPtrIn,
36  ParticleData* particleDataPtrIn, Couplings* couplingsPtrIn) {
37 
38  // Save pointers.
39  infoPtr = infoPtrIn;
40  settingsPtr = settingsPtrIn;
41  particleDataPtr = particleDataPtrIn;
42  couplingsPtr = couplingsPtrIn;
43 
44  // Minimal decaying-resonance width. Minimal phase space for meMode = 103.
45  minWidth = settingsPtr->parm("ResonanceWidths:minWidth");
46  minThreshold = settingsPtr->parm("ResonanceWidths:minThreshold");
47 
48  // Pointer to particle species.
49  particlePtr = particleDataPtr->particleDataEntryPtr(idRes);
50  if (particlePtr == 0) {
51  infoPtr->errorMsg("Error in ResonanceWidths::init:"
52  " unknown resonance identity code");
53  return false;
54  }
55 
56  // Generic particles should not have meMode < 100, but allow
57  // some exceptions: not used Higgses and not used Technicolor.
58  if (idRes == 35 || idRes == 36 || idRes == 37
59  || idRes/1000000 == 3) isGeneric = false;
60 
61  // Resonance properties: antiparticle, mass, width
62  hasAntiRes = particlePtr->hasAnti();
63  mRes = particlePtr->m0();
64  GammaRes = particlePtr->mWidth();
65  m2Res = mRes*mRes;
66 
67  // For very narrow resonances assign fictitious small width.
68  if (GammaRes < minWidth) GammaRes = 0.1 * minWidth;
69  GamMRat = GammaRes / mRes;
70 
71  // Secondary decay chains by default all on.
72  openPos = 1.;
73  openNeg = 1.;
74 
75  // Allow option where on-shell width is forced to current value.
76  doForceWidth = particlePtr->doForceWidth();
77  forceFactor = 1.;
78 
79  // Check that resonance OK.
80  if (particlePtr == 0) infoPtr->errorMsg("Error in ResonanceWidths::init:"
81  " unknown resonance identity code");
82 
83  // Initialize constants used for a resonance.
84  initConstants();
85 
86  // Calculate various common prefactors for the current mass.
87  mHat = mRes;
88  calcPreFac(true);
89 
90  // Reset quantities to sum. Declare variables inside loop.
91  double widTot = 0.;
92  double widPos = 0.;
93  double widNeg = 0.;
94  int idNow, idAnti;
95  double openSecPos, openSecNeg;
96 
97  // Loop over all decay channels. Basic properties of channel.
98  for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
99  iChannel = i;
100  onMode = particlePtr->channel(i).onMode();
101  meMode = particlePtr->channel(i).meMode();
102  mult = particlePtr->channel(i).multiplicity();
103  widNow = 0.;
104 
105  // Warn if not relevant meMode.
106  if ( meMode < 0 || meMode > 103 || (isGeneric && meMode < 100) ) {
107  infoPtr->errorMsg("Error in ResonanceWidths::init:"
108  " resonance meMode not acceptable");
109  }
110 
111  // Channels with meMode < 100 must be implemented in derived classes.
112  if (meMode < 100) {
113 
114  // Read out information on channel: primarily use first two.
115  id1 = particlePtr->channel(i).product(0);
116  id2 = particlePtr->channel(i).product(1);
117  id1Abs = abs(id1);
118  id2Abs = abs(id2);
119 
120  // Order first two in descending order of absolute values.
121  if (id2Abs > id1Abs) {swap( id1, id2); swap( id1Abs, id2Abs);}
122 
123  // Allow for third product to be treated in derived classes.
124  if (mult > 2) {
125  id3 = particlePtr->channel(i).product(2);
126  id3Abs = abs(id3);
127 
128  // Also order third into descending order of absolute values.
129  if (id3Abs > id2Abs) {swap( id2, id3); swap( id2Abs, id3Abs);}
130  if (id2Abs > id1Abs) {swap( id1, id2); swap( id1Abs, id2Abs);}
131  }
132 
133  // Read out masses. Calculate two-body phase space.
134  mf1 = particleDataPtr->m0(id1Abs);
135  mf2 = particleDataPtr->m0(id2Abs);
136  mr1 = pow2(mf1 / mHat);
137  mr2 = pow2(mf2 / mHat);
138  ps = (mHat < mf1 + mf2 + MASSMARGIN) ? 0.
139  : sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 );
140  if (mult > 2) {
141  mf3 = particleDataPtr->m0(id3Abs);
142  mr3 = pow2(mf3 / mHat);
143  }
144 
145  // Let derived class calculate width for channel provided.
146  calcWidth(true);
147  }
148 
149  // Channels with meMode >= 100 are calculated based on stored values.
150  else widNow = GammaRes * particlePtr->channel(i).bRatio();
151 
152  // Find secondary open fractions of partial width.
153  openSecPos = 1.;
154  openSecNeg = 1.;
155  if (widNow > 0.) for (int j = 0; j < mult; ++j) {
156  idNow = particlePtr->channel(i).product(j);
157  idAnti = (particleDataPtr->hasAnti(idNow)) ? -idNow : idNow;
158  // Secondary widths not yet initialized for heavier states,
159  // so have to assume unit open fraction there.
160  if (idNow == 23 || abs(idNow) == 24
161  || particleDataPtr->m0(abs(idNow)) < mRes) {
162  openSecPos *= particleDataPtr->resOpenFrac(idNow);
163  openSecNeg *= particleDataPtr->resOpenFrac(idAnti);
164  }
165  }
166 
167  // Store partial widths and secondary open fractions.
168  particlePtr->channel(i).onShellWidth(widNow);
169  particlePtr->channel(i).openSec( idRes, openSecPos);
170  particlePtr->channel(i).openSec(-idRes, openSecNeg);
171 
172  // Update sum over all channnels and over open channels only.
173  widTot += widNow;
174  if (onMode == 1 || onMode == 2) widPos += widNow * openSecPos;
175  if (onMode == 1 || onMode == 3) widNeg += widNow * openSecNeg;
176  }
177 
178  // If no decay channels are open then set particle stable and done.
179  if (widTot < minWidth) {
180  particlePtr->setMayDecay(false, false);
181  particlePtr->setMWidth(0., false);
182  for (int i = 0; i < particlePtr->sizeChannels(); ++i)
183  particlePtr->channel(i).bRatio( 0., false);
184  return true;
185  }
186 
187  // Normalize branching ratios to unity.
188  double bRatio;
189  for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
190  bRatio = particlePtr->channel(i).onShellWidth() / widTot;
191  particlePtr->channel(i).bRatio( bRatio, false);
192  }
193 
194  // Optionally force total width by rescaling of all partial ones.
195  if (doForceWidth) {
196  forceFactor = GammaRes / widTot;
197  for (int i = 0; i < particlePtr->sizeChannels(); ++i)
198  particlePtr->channel(i).onShellWidthFactor( forceFactor);
199  }
200 
201  // Else update newly calculated partial width.
202  else {
203  particlePtr->setMWidth(widTot, false);
204  GammaRes = widTot;
205  }
206 
207  // Updated width-to-mass ratio. Secondary widths for open.
208  GamMRat = GammaRes / mRes;
209  openPos = widPos / widTot;
210  openNeg = widNeg / widTot;
211 
212  // Clip wings of Higgses.
213  bool isHiggs = (idRes == 25 || idRes == 35 ||idRes == 36 ||idRes == 37);
214  bool clipHiggsWings = settingsPtr->flag("Higgs:clipWings");
215  if (isHiggs && clipHiggsWings) {
216  double mMinNow = particlePtr->mMin();
217  double mMaxNow = particlePtr->mMax();
218  double wingsFac = settingsPtr->parm("Higgs:wingsFac");
219  double mMinWing = mRes - wingsFac * GammaRes;
220  double mMaxWing = mRes + wingsFac * GammaRes;
221  if (mMinWing > mMinNow) particlePtr->setMMinNoChange(mMinWing);
222  if (mMaxWing < mMaxNow || mMaxNow < mMinNow)
223  particlePtr->setMMaxNoChange(mMaxWing);
224  }
225 
226  // Done.
227  return true;
228 
229 }
230 
231 //--------------------------------------------------------------------------
232 
233 // Calculate the total width and store phase-space-weighted coupling sums.
234 
235 double ResonanceWidths::width(int idSgn, double mHatIn, int idInFlavIn,
236  bool openOnly, bool setBR, int idOutFlav1, int idOutFlav2) {
237 
238  // Calculate various prefactors for the current mass.
239  mHat = mHatIn;
240  idInFlav = idInFlavIn;
241  calcPreFac(false);
242 
243  // Reset quantities to sum. Declare variables inside loop.
244  double widSum = 0.;
245  double mfSum, psOnShell;
246 
247  // Loop over all decay channels. Basic properties of channel.
248  for (int i = 0; i < particlePtr->sizeChannels(); ++i) {
249  iChannel = i;
250  onMode = particlePtr->channel(i).onMode();
251  meMode = particlePtr->channel(i).meMode();
252  mult = particlePtr->channel(i).multiplicity();
253 
254  // Initially assume vanishing branching ratio.
255  widNow = 0.;
256  if (setBR) particlePtr->channel(i).currentBR(widNow);
257 
258  // Optionally only consider specific (two-body) decay channel.
259  // Currently only used for Higgs -> q qbar, g g or gamma gamma.
260  if (idOutFlav1 > 0 || idOutFlav2 > 0) {
261  if (mult > 2) continue;
262  if (particlePtr->channel(i).product(0) != idOutFlav1) continue;
263  if (particlePtr->channel(i).product(1) != idOutFlav2) continue;
264  }
265 
266  // Optionally only consider open channels.
267  if (openOnly) {
268  if (idSgn > 0 && onMode !=1 && onMode != 2) continue;
269  if (idSgn < 0 && onMode !=1 && onMode != 3) continue;
270  }
271 
272  // Channels with meMode < 100 must be implemented in derived classes.
273  if (meMode < 100) {
274 
275  // Read out information on channel: primarily use first two.
276  id1 = particlePtr->channel(i).product(0);
277  id2 = particlePtr->channel(i).product(1);
278  id1Abs = abs(id1);
279  id2Abs = abs(id2);
280 
281  // Order first two in descending order of absolute values.
282  if (id2Abs > id1Abs) {swap( id1, id2); swap( id1Abs, id2Abs);}
283 
284  // Allow for third product to be treated in derived classes.
285  if (mult > 2) {
286  id3 = particlePtr->channel(i).product(2);
287  id3Abs = abs(id3);
288 
289  // Also order third into descending order of absolute values.
290  if (id3Abs > id2Abs) {swap( id2, id3); swap( id2Abs, id3Abs);}
291  if (id2Abs > id1Abs) {swap( id1, id2); swap( id1Abs, id2Abs);}
292  }
293 
294  // Read out masses. Calculate two-body phase space.
295  mf1 = particleDataPtr->m0(id1Abs);
296  mf2 = particleDataPtr->m0(id2Abs);
297  mr1 = pow2(mf1 / mHat);
298  mr2 = pow2(mf2 / mHat);
299  ps = (mHat < mf1 + mf2 + MASSMARGIN) ? 0.
300  : sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 );
301  if (mult > 2) {
302  mf3 = particleDataPtr->m0(id3Abs);
303  mr3 = pow2(mf3 / mHat);
304  }
305 
306  // Let derived class calculate width for channel provided.
307  calcWidth(false);
308  }
309 
310  // Now on to meMode >= 100. First case: no correction at all.
311  else if (meMode == 100)
312  widNow = GammaRes * particlePtr->channel(i).bRatio();
313 
314  // Correction by step at threshold.
315  else if (meMode == 101) {
316  mfSum = 0.;
317  for (int j = 0; j < mult; ++j) mfSum
318  += particleDataPtr->m0( particlePtr->channel(i).product(j) );
319  if (mfSum + MASSMARGIN < mHat)
320  widNow = GammaRes * particlePtr->channel(i).bRatio();
321  }
322 
323  // Correction by a phase space factor for two-body decays.
324  else if ( (meMode == 102 || meMode == 103) && mult == 2) {
325  mf1 = particleDataPtr->m0( particlePtr->channel(i).product(0) );
326  mf2 = particleDataPtr->m0( particlePtr->channel(i).product(1) );
327  mr1 = pow2(mf1 / mHat);
328  mr2 = pow2(mf2 / mHat);
329  ps = (mHat < mf1 + mf2 + MASSMARGIN) ? 0.
330  : sqrtpos( pow2(1. - mr1 - mr2) - 4. * mr1 * mr2 );
331  mr1 = pow2(mf1 / mRes);
332  mr2 = pow2(mf2 / mRes);
333  psOnShell = (meMode == 102) ? 1. : max( minThreshold,
334  sqrtpos( pow2(1.- mr1 - mr2) - 4. * mr1 * mr2) );
335  widNow = GammaRes * particlePtr->channel(i).bRatio() * ps / psOnShell;
336  }
337 
338  // Correction by simple threshold factor for multibody decay.
339  else if (meMode == 102 || meMode == 103) {
340  mfSum = 0.;
341  for (int j = 0; j < mult; ++j) mfSum
342  += particleDataPtr->m0( particlePtr->channel(i).product(j) );
343  ps = sqrtpos(1. - mfSum / mHat);
344  psOnShell = (meMode == 102) ? 1. : max( minThreshold,
345  sqrtpos(1. - mfSum / mRes) );
346  widNow = GammaRes * particlePtr->channel(i).bRatio() * ps / psOnShell;
347  }
348 
349  // Optionally multiply by secondary widths.
350  if (openOnly) widNow *= particlePtr->channel(i).openSec(idSgn);
351 
352  // Optionally include factor to force to fixed width??
353  // Optionally multiply by current/nominal resonance mass??
354 
355  // Sum back up.
356  widSum += widNow;
357 
358  // Optionally store partial widths for later decay channel choice.
359  if (setBR) particlePtr->channel(i).currentBR(widNow);
360  }
361 
362  // Done.
363  return widSum;
364 
365 }
366 
367 //--------------------------------------------------------------------------
368 
369 // Numerical integration of matrix-element in two-body decay,
370 // where one particle is described by a Breit-Wigner mass distribution.
371 // Normalization to unit integral if matrix element is unity
372 // and there are no phase-space restrictions.
373 
374 double ResonanceWidths::numInt1BW(double mHatIn, double m1, double Gamma1,
375  double mMin1, double m2, int psMode) {
376 
377  // Check that phase space is open for integration.
378  if (mMin1 + m2 > mHatIn) return 0.;
379 
380  // Precalculate coefficients for Breit-Wigner selection.
381  double s1 = m1 * m1;
382  double mG1 = m1 * Gamma1;
383  double mMax1 = mHatIn - m2;
384  double atanMin1 = atan( (mMin1 * mMin1 - s1) / mG1 );
385  double atanMax1 = atan( (mMax1 * mMax1 - s1) / mG1 );
386  double atanDif1 = atanMax1 - atanMin1;
387  double wtDif1 = atanDif1 / (M_PI * NPOINT);
388 
389  // Step size in atan-mapped variable.
390  double xStep = 1. / NPOINT;
391 
392  // Variables used in loop over integration points.
393  double sum = 0.;
394  double mrNow2 = pow2(m2 / mHatIn);
395  double xNow1, sNow1, mNow1, mrNow1, psNow, value;
396 
397  // Loop with first-particle mass selection.
398  for (int ip1 = 0; ip1 < NPOINT; ++ip1) {
399  xNow1 = xStep * (ip1 + 0.5);
400  sNow1 = s1 + mG1 * tan(atanMin1 + xNow1 * atanDif1);
401  mNow1 = min( mMax1, max( mMin1, sqrtpos(sNow1) ) );
402  mrNow1 = pow2(mNow1 / mHatIn);
403 
404  // Evaluate value and add to sum. Different matrix elements.
405  psNow = sqrtpos( pow2(1. - mrNow1 - mrNow2)
406  - 4. * mrNow1 * mrNow2);
407  value = 1.;
408  if (psMode == 1) value = psNow;
409  else if (psMode == 2) value = psNow * psNow;
410  else if (psMode == 3) value = pow3(psNow);
411  else if (psMode == 5) value = psNow
412  * (pow2(1. - mrNow1 - mrNow2) + 8. * mrNow1 * mrNow2);
413  else if (psMode == 6) value = pow3(psNow);
414  sum += value;
415 
416  // End of loop over integration points. Overall normalization.
417  }
418  sum *= wtDif1;
419 
420  // Done.
421  return sum;
422 }
423 
424 //--------------------------------------------------------------------------
425 
426 // Numerical integration of matrix-element in two-body decay,
427 // where both particles are described by Breit-Wigner mass distributions.
428 // Normalization to unit integral if matrix element is unity
429 // and there are no phase-space restrictions.
430 
431 double ResonanceWidths::numInt2BW(double mHatIn, double m1, double Gamma1,
432  double mMin1, double m2, double Gamma2, double mMin2, int psMode) {
433 
434  // Check that phase space is open for integration.
435  if (mMin1 + mMin2 >= mHatIn) return 0.;
436 
437  // Precalculate coefficients for Breit-Wigner selection.
438  double s1 = m1 * m1;
439  double mG1 = m1 * Gamma1;
440  double mMax1 = mHatIn - mMin2;
441  double atanMin1 = atan( (mMin1 * mMin1 - s1) / mG1 );
442  double atanMax1 = atan( (mMax1 * mMax1 - s1) / mG1 );
443  double atanDif1 = atanMax1 - atanMin1;
444  double wtDif1 = atanDif1 / (M_PI * NPOINT);
445  double s2 = m2 * m2;
446  double mG2 = m2 * Gamma2;
447  double mMax2 = mHatIn - mMin1;
448  double atanMin2 = atan( (mMin2 * mMin2 - s2) / mG2 );
449  double atanMax2 = atan( (mMax2 * mMax2 - s2) / mG2 );
450  double atanDif2 = atanMax2 - atanMin2;
451  double wtDif2 = atanDif2 / (M_PI * NPOINT);
452 
453  // If on-shell decay forbidden then split integration range
454  // to ensure that low-mass region is not forgotten.
455  bool mustDiv = false;
456  double mDiv1 = 0.;
457  double atanDiv1 = 0.;
458  double atanDLo1 = 0.;
459  double atanDHi1 = 0.;
460  double wtDLo1 = 0.;
461  double wtDHi1 = 0.;
462  double mDiv2 = 0.;
463  double atanDiv2 = 0.;
464  double atanDLo2 = 0.;
465  double atanDHi2 = 0.;
466  double wtDLo2 = 0.;
467  double wtDHi2 = 0.;
468  if (m1 + m2 > mHatIn) {
469  mustDiv = true;
470  double tmpDiv = (mHatIn - m1 - m2) / (Gamma1 + Gamma2);
471  mDiv1 = m1 + Gamma1 * tmpDiv;
472  atanDiv1 = atan( (mDiv1 * mDiv1 - s1) / mG1 );
473  atanDLo1 = atanDiv1 - atanMin1;
474  atanDHi1 = atanMax1 - atanDiv1;
475  wtDLo1 = atanDLo1 / (M_PI * NPOINT);
476  wtDHi1 = atanDHi1 / (M_PI * NPOINT);
477  mDiv2 = m2 + Gamma2 * tmpDiv;
478  atanDiv2 = atan( (mDiv2 * mDiv2 - s2) / mG2 );
479  atanDLo2 = atanDiv2 - atanMin2;
480  atanDHi2 = atanMax2 - atanDiv2;
481  wtDLo2 = atanDLo2 / (M_PI * NPOINT);
482  wtDHi2 = atanDHi2 / (M_PI * NPOINT);
483  }
484 
485  // Step size in atan-mapped variable.
486  double xStep = 1. / NPOINT;
487  int nIter = (mustDiv) ? 2 * NPOINT : NPOINT;
488 
489  // Variables used in loop over integration points.
490  double sum = 0.;
491  double xNow1, sNow1, mNow1, mrNow1, xNow2, sNow2, mNow2, mrNow2, psNow,
492  value;
493  double wtNow1 = wtDif1;
494  double wtNow2 = wtDif2;
495 
496  // Outer loop with first-particle mass selection.
497  for (int ip1 = 0; ip1 < nIter; ++ip1) {
498  if (!mustDiv) {
499  xNow1 = xStep * (ip1 + 0.5);
500  sNow1 = s1 + mG1 * tan(atanMin1 + xNow1 * atanDif1);
501  } else if (ip1 < NPOINT) {
502  xNow1 = xStep * (ip1 + 0.5);
503  sNow1 = s1 + mG1 * tan(atanMin1 + xNow1 * atanDLo1);
504  wtNow1 = wtDLo1;
505  } else {
506  xNow1 = xStep * (ip1 - NPOINT + 0.5);
507  sNow1 = s1 + mG1 * tan(atanDiv1 + xNow1 * atanDHi1);
508  wtNow1 = wtDHi1;
509  }
510  mNow1 = min( mMax1, max( mMin1, sqrtpos(sNow1) ) );
511  mrNow1 = pow2(mNow1 / mHatIn);
512 
513  // Inner loop with second-particle mass selection.
514  for (int ip2 = 0; ip2 < nIter; ++ip2) {
515  if (!mustDiv) {
516  xNow2 = xStep * (ip2 + 0.5);
517  sNow2 = s2 + mG2 * tan(atanMin2 + xNow2 * atanDif2);
518  } else if (ip2 < NPOINT) {
519  xNow2 = xStep * (ip2 + 0.5);
520  sNow2 = s2 + mG2 * tan(atanMin2 + xNow2 * atanDLo2);
521  wtNow2 = wtDLo2;
522  } else {
523  xNow2 = xStep * (ip2 - NPOINT + 0.5);
524  sNow2 = s2 + mG2 * tan(atanDiv2 + xNow2 * atanDHi2);
525  wtNow2 = wtDHi2;
526  }
527  mNow2 = min( mMax2, max( mMin2, sqrtpos(sNow2) ) );
528  mrNow2 = pow2(mNow2 / mHatIn);
529 
530  // Check that point is inside phase space.
531  if (mNow1 + mNow2 > mHatIn) break;
532 
533  // Evaluate value and add to sum. Different matrix elements.
534  psNow = sqrtpos( pow2(1. - mrNow1 - mrNow2)
535  - 4. * mrNow1 * mrNow2);
536  value = 1.;
537  if (psMode == 1) value = psNow;
538  else if (psMode == 2) value = psNow * psNow;
539  else if (psMode == 3) value = pow3(psNow);
540  else if (psMode == 5) value = psNow
541  * (pow2(1. - mrNow1 - mrNow2) + 8. * mrNow1 * mrNow2);
542  else if (psMode == 6) value = pow3(psNow);
543  sum += value * wtNow1 * wtNow2;
544 
545  // End of second and first loop over integration points.
546  }
547  }
548 
549  // Done.
550  return sum;
551 }
552 
553 //==========================================================================
554 
555 // The ResonanceGmZ class.
556 // Derived class for gamma*/Z0 properties.
557 
558 //--------------------------------------------------------------------------
559 
560 // Initialize constants.
561 
562 void ResonanceGmZ::initConstants() {
563 
564  // Locally stored properties and couplings.
565  gmZmode = settingsPtr->mode("WeakZ0:gmZmode");
566  thetaWRat = 1. / (16. * couplingsPtr->sin2thetaW()
567  * couplingsPtr->cos2thetaW());
568 
569 }
570 
571 //--------------------------------------------------------------------------
572 
573 // Calculate various common prefactors for the current mass.
574 
575 void ResonanceGmZ::calcPreFac(bool calledFromInit) {
576 
577  // Common coupling factors.
578  alpEM = couplingsPtr->alphaEM(mHat * mHat);
579  alpS = couplingsPtr->alphaS(mHat * mHat);
580  colQ = 3. * (1. + alpS / M_PI);
581  preFac = alpEM * thetaWRat * mHat / 3.;
582 
583  // When call for incoming flavour need to consider gamma*/Z0 mix.
584  if (!calledFromInit) {
585 
586  // Couplings when an incoming fermion is specified; elso only pure Z0.
587  ei2 = 0.;
588  eivi = 0.;
589  vi2ai2 = 1.;
590  int idInFlavAbs = abs(idInFlav);
591  if (idInFlavAbs > 0 && idInFlavAbs < 19) {
592  ei2 = couplingsPtr->ef2(idInFlavAbs);
593  eivi = couplingsPtr->efvf(idInFlavAbs);
594  vi2ai2 = couplingsPtr->vf2af2(idInFlavAbs);
595  }
596 
597  // Calculate prefactors for gamma/interference/Z0 terms.
598  double sH = mHat * mHat;
599  gamNorm = ei2;
600  intNorm = 2. * eivi * thetaWRat * sH * (sH - m2Res)
601  / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
602  resNorm = vi2ai2 * pow2(thetaWRat * sH)
603  / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
604 
605  // Rescale Z0 height normalization to compensate for a width one??
606  //if (doForceWidth) {
607  // intNorm *= forceFactor;
608  // resNorm *= forceFactor;
609  //}
610 
611  // Optionally only keep gamma* or Z0 term.
612  if (gmZmode == 1) {intNorm = 0.; resNorm = 0.;}
613  if (gmZmode == 2) {gamNorm = 0.; intNorm = 0.;}
614  }
615 
616 }
617 
618 //--------------------------------------------------------------------------
619 
620 // Calculate width for currently considered channel.
621 
622 void ResonanceGmZ::calcWidth(bool calledFromInit) {
623 
624  // Check that above threshold.
625  if (ps == 0.) return;
626 
627  // Only contributions from three fermion generations, except top.
628  if ( (id1Abs > 5 && id1Abs < 11) || id1Abs > 16 ) return;
629 
630  // At initialization only the pure Z0 should be considered.
631  if (calledFromInit) {
632 
633  // Combine kinematics with colour factor and couplings.
634  widNow = preFac * ps * (couplingsPtr->vf2(id1Abs) * (1. + 2. * mr1)
635  + couplingsPtr->af2(id1Abs) * ps*ps);
636  if (id1Abs < 6) widNow *= colQ;
637  }
638 
639  // When call for incoming flavour need to consider gamma*/Z0 mix.
640  else {
641 
642  // Kinematical factors and couplings.
643  double kinFacV = ps * (1. + 2. * mr1);
644  double ef2 = couplingsPtr->ef2(id1Abs) * kinFacV;
645  double efvf = couplingsPtr->efvf(id1Abs) * kinFacV;
646  double vf2af2 = couplingsPtr->vf2(id1Abs) * kinFacV
647  + couplingsPtr->af2(id1Abs) * pow3(ps);
648 
649  // Relative outwidths: combine instate, propagator and outstate.
650  widNow = gamNorm * ef2 + intNorm * efvf + resNorm * vf2af2;
651 
652  // Colour factor.
653  if (id1Abs < 6) widNow *= colQ;
654  }
655 
656 }
657 
658 //==========================================================================
659 
660 // The ResonanceW class.
661 // Derived class for W+- properties.
662 
663 //--------------------------------------------------------------------------
664 
665 // Initialize constants.
666 
667 void ResonanceW::initConstants() {
668 
669  // Locally stored properties and couplings.
670  thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());
671 
672 }
673 
674 //--------------------------------------------------------------------------
675 
676 // Calculate various common prefactors for the current mass.
677 
678 void ResonanceW::calcPreFac(bool) {
679 
680  // Common coupling factors.
681  alpEM = couplingsPtr->alphaEM(mHat * mHat);
682  alpS = couplingsPtr->alphaS(mHat * mHat);
683  colQ = 3. * (1. + alpS / M_PI);
684  preFac = alpEM * thetaWRat * mHat;
685 
686 }
687 
688 //--------------------------------------------------------------------------
689 
690 // Calculate width for currently considered channel.
691 
692 void ResonanceW::calcWidth(bool) {
693 
694  // Check that above threshold.
695  if (ps == 0.) return;
696 
697  // Only contributions from three fermion generations, except top.
698  if ( (id1Abs > 5 && id1Abs < 11) || id1Abs > 16 ) return;
699 
700 
701  // Combine kinematics with colour factor and couplings.
702  widNow = preFac * ps
703  * (1. - 0.5 * (mr1 + mr2) - 0.5 * pow2(mr1 - mr2));
704  if (id1Abs < 6) widNow *= colQ * couplingsPtr->V2CKMid(id1Abs, id2Abs);
705 
706 }
707 
708 //==========================================================================
709 
710 // The ResonanceTop class.
711 // Derived class for top/antitop properties.
712 
713 //--------------------------------------------------------------------------
714 
715 // Initialize constants.
716 
717 void ResonanceTop::initConstants() {
718 
719  // Locally stored properties and couplings.
720  thetaWRat = 1. / (16. * couplingsPtr->sin2thetaW());
721  m2W = pow2(particleDataPtr->m0(24));
722 
723 }
724 
725 //--------------------------------------------------------------------------
726 
727 // Calculate various common prefactors for the current mass.
728 
729 void ResonanceTop::calcPreFac(bool) {
730 
731  // Common coupling factors.
732  alpEM = couplingsPtr->alphaEM(mHat * mHat);
733  alpS = couplingsPtr->alphaS(mHat * mHat);
734  colQ = 1. - 2.5 * alpS / M_PI;
735  preFac = alpEM * thetaWRat * pow3(mHat) / m2W;
736 
737 }
738 
739 //--------------------------------------------------------------------------
740 
741 // Calculate width for currently considered channel.
742 
743 void ResonanceTop::calcWidth(bool) {
744 
745  // Only contributions from W + quark.
746  if (id1Abs != 24 || id2Abs > 5) return;
747 
748  // Check that above threshold. Kinematical factor.
749  if (ps == 0.) return;
750  widNow = preFac * ps
751  * ( pow2(1. - mr2) + (1. + mr2) * mr1 - 2. * mr1 * mr1 );
752 
753  // Combine with colour factor and CKM couplings.
754  widNow *= colQ * couplingsPtr->V2CKMid(6, id2Abs);
755 
756 }
757 
758 //==========================================================================
759 
760 // The ResonanceFour class.
761 // Derived class for fourth-generation properties.
762 
763 //--------------------------------------------------------------------------
764 
765 // Initialize constants.
766 
767 void ResonanceFour::initConstants() {
768 
769  // Locally stored properties and couplings.
770  thetaWRat = 1. / (16. * couplingsPtr->sin2thetaW());
771  m2W = pow2(particleDataPtr->m0(24));
772 
773 }
774 
775 //--------------------------------------------------------------------------
776 
777 // Calculate various common prefactors for the current mass.
778 
779 void ResonanceFour::calcPreFac(bool) {
780 
781  // Common coupling factors.
782  alpEM = couplingsPtr->alphaEM(mHat * mHat);
783  alpS = couplingsPtr->alphaS(mHat * mHat);
784  colQ = (idRes < 9) ? 1. - 2.5 * alpS / M_PI : 1.;
785  preFac = alpEM * thetaWRat * pow3(mHat) / m2W;
786 
787 }
788 
789 //--------------------------------------------------------------------------
790 
791 // Calculate width for currently considered channel.
792 
793 void ResonanceFour::calcWidth(bool) {
794 
795  // Only contributions from W + fermion.
796  if (id1Abs != 24 || id2Abs > 18) return;
797 
798  // Check that above threshold. Kinematical factor.
799  if (ps == 0.) return;
800  widNow = preFac * ps
801  * ( pow2(1. - mr2) + (1. + mr2) * mr1 - 2. * mr1 * mr1 );
802 
803  // Combine with colour factor and CKM couplings.
804  if (idRes < 9) widNow *= colQ * couplingsPtr->V2CKMid(idRes, id2Abs);
805 
806 }
807 
808 //==========================================================================
809 
810 // The ResonanceH class.
811 // Derived class for SM and BSM Higgs properties.
812 
813 //--------------------------------------------------------------------------
814 
815 // Constants: could be changed here if desired, but normally should not.
816 // These are of technical nature, as described for each.
817 
818 // Minimal mass for W, Z, top in integration over respective Breit-Wigner.
819 // Top constrainted by t -> W b decay, which is not seen in simple top BW.
820 const double ResonanceH::MASSMINWZ = 10.;
821 const double ResonanceH::MASSMINT = 100.;
822 
823 // Number of widths above threshold where B-W integration not needed.
824 const double ResonanceH::GAMMAMARGIN = 10.;
825 
826 //--------------------------------------------------------------------------
827 
828 // Initialize constants.
829 
830 void ResonanceH::initConstants() {
831 
832  // Locally stored properties and couplings.
833  useCubicWidth = settingsPtr->flag("Higgs:cubicWidth");
834  useRunLoopMass = settingsPtr->flag("Higgs:runningLoopMass");
835  sin2tW = couplingsPtr->sin2thetaW();
836  cos2tW = 1. - sin2tW;
837  mT = particleDataPtr->m0(6);
838  mZ = particleDataPtr->m0(23);
839  mW = particleDataPtr->m0(24);
840  mHchg = particleDataPtr->m0(37);
841  GammaT = particleDataPtr->mWidth(6);
842  GammaZ = particleDataPtr->mWidth(23);
843  GammaW = particleDataPtr->mWidth(24);
844 
845  // Couplings to fermions, Z and W, depending on Higgs type.
846  coup2d = 1.;
847  coup2u = 1.;
848  coup2l = 1.;
849  coup2Z = 1.;
850  coup2W = 1.;
851  coup2Hchg = 0.;
852  coup2H1H1 = 0.;
853  coup2A3A3 = 0.;
854  coup2H1Z = 0.;
855  coup2A3Z = 0.;
856  coup2A3H1 = 0.;
857  coup2HchgW = 0.;
858  if (higgsType == 1) {
859  coup2d = settingsPtr->parm("HiggsH1:coup2d");
860  coup2u = settingsPtr->parm("HiggsH1:coup2u");
861  coup2l = settingsPtr->parm("HiggsH1:coup2l");
862  coup2Z = settingsPtr->parm("HiggsH1:coup2Z");
863  coup2W = settingsPtr->parm("HiggsH1:coup2W");
864  coup2Hchg = settingsPtr->parm("HiggsH1:coup2Hchg");
865  } else if (higgsType == 2) {
866  coup2d = settingsPtr->parm("HiggsH2:coup2d");
867  coup2u = settingsPtr->parm("HiggsH2:coup2u");
868  coup2l = settingsPtr->parm("HiggsH2:coup2l");
869  coup2Z = settingsPtr->parm("HiggsH2:coup2Z");
870  coup2W = settingsPtr->parm("HiggsH2:coup2W");
871  coup2Hchg = settingsPtr->parm("HiggsH2:coup2Hchg");
872  coup2H1H1 = settingsPtr->parm("HiggsH2:coup2H1H1");
873  coup2A3A3 = settingsPtr->parm("HiggsH2:coup2A3A3");
874  coup2H1Z = settingsPtr->parm("HiggsH2:coup2H1Z");
875  coup2A3Z = settingsPtr->parm("HiggsA3:coup2H2Z");
876  coup2A3H1 = settingsPtr->parm("HiggsH2:coup2A3H1");
877  coup2HchgW = settingsPtr->parm("HiggsH2:coup2HchgW");
878  } else if (higgsType == 3) {
879  coup2d = settingsPtr->parm("HiggsA3:coup2d");
880  coup2u = settingsPtr->parm("HiggsA3:coup2u");
881  coup2l = settingsPtr->parm("HiggsA3:coup2l");
882  coup2Z = settingsPtr->parm("HiggsA3:coup2Z");
883  coup2W = settingsPtr->parm("HiggsA3:coup2W");
884  coup2Hchg = settingsPtr->parm("HiggsA3:coup2Hchg");
885  coup2H1H1 = settingsPtr->parm("HiggsA3:coup2H1H1");
886  coup2H1Z = settingsPtr->parm("HiggsA3:coup2H1Z");
887  coup2HchgW = settingsPtr->parm("HiggsA3:coup2Hchg");
888  }
889 
890  // Initialization of threshold kinematical factor by stepwise
891  // numerical integration of H -> t tbar, Z0 Z0 and W+ W-.
892  int psModeT = (higgsType < 3) ? 3 : 1;
893  int psModeWZ = (higgsType < 3) ? 5 : 6;
894  mLowT = max( 2.02 * MASSMINT, 0.5 * mT);
895  mStepT = 0.01 * (3. * mT - mLowT);
896  mLowZ = max( 2.02 * MASSMINWZ, 0.5 * mZ);
897  mStepZ = 0.01 * (3. * mZ - mLowZ);
898  mLowW = max( 2.02 * MASSMINWZ, 0.5 * mW);
899  mStepW = 0.01 * (3. * mW - mLowW);
900  for (int i = 0; i <= 100; ++i) {
901  kinFacT[i] = numInt2BW( mLowT + i * mStepT,
902  mT, GammaT, MASSMINT, mT, GammaT, MASSMINT, psModeT);
903  kinFacZ[i] = numInt2BW( mLowZ + i * mStepZ,
904  mZ, GammaZ, MASSMINWZ, mZ, GammaZ, MASSMINWZ, psModeWZ);
905  kinFacW[i] = numInt2BW( mLowW + i * mStepW,
906  mW, GammaW, MASSMINWZ, mW, GammaW, MASSMINWZ, psModeWZ);
907  }
908 
909 }
910 
911 //--------------------------------------------------------------------------
912 
913 // Calculate various common prefactors for the current mass.
914 
915 void ResonanceH::calcPreFac(bool) {
916 
917  // Common coupling factors.
918  alpEM = couplingsPtr->alphaEM(mHat * mHat);
919  alpS = couplingsPtr->alphaS(mHat * mHat);
920  colQ = 3. * (1. + alpS / M_PI);
921  preFac = (alpEM / (8. * sin2tW)) * pow3(mHat) / pow2(mW);
922 }
923 
924 //--------------------------------------------------------------------------
925 
926 // Calculate width for currently considered channel.
927 
928 void ResonanceH::calcWidth(bool) {
929 
930  // Widths of decays Higgs -> f + fbar.
931  if ( id2Abs == id1Abs && ( (id1Abs > 0 && id1Abs < 7)
932  || (id1Abs > 10 && id1Abs < 17) ) ) {
933  kinFac = 0.;
934 
935  // Check that above threshold (well above for top). Kinematical factor.
936  if ( (id1Abs != 6 && mHat > 2. * mf1 + MASSMARGIN)
937  || (id1Abs == 6 && mHat > 3. * mT ) ) {
938  // A0 behaves like beta, h0 and H0 like beta**3.
939  kinFac = (higgsType < 3) ? pow3(ps) : ps;
940  }
941 
942  // Top near or below threshold: interpolate in table.
943  else if (id1Abs == 6 && mHat > mLowT) {
944  double xTab = (mHat - mLowT) / mStepT;
945  int iTab = max( 0, min( 99, int(xTab) ) );
946  kinFac = kinFacT[iTab]
947  * pow( kinFacT[iTab + 1] / kinFacT[iTab], xTab - iTab);
948  }
949 
950  // Coupling from mass and from BSM deviation from SM.
951  double coupFac = pow2(particleDataPtr->mRun(id1Abs, mHat) / mHat);
952  if (id1Abs < 7 && id1Abs%2 == 1) coupFac *= coup2d * coup2d;
953  else if (id1Abs < 7) coupFac *= coup2u * coup2u;
954  else coupFac *= coup2l * coup2l;
955 
956  // Combine couplings and phase space with colour factor.
957  widNow = preFac * coupFac * kinFac;
958  if (id1Abs < 7) widNow *= colQ;
959  }
960 
961  // Widths of decays Higgs -> g + g.
962  else if (id1Abs == 21 && id2Abs == 21)
963  widNow = preFac * pow2(alpS / M_PI) * eta2gg();
964 
965  // Widths of decays Higgs -> gamma + gamma.
966  else if (id1Abs == 22 && id2Abs == 22)
967  widNow = preFac * pow2(alpEM / M_PI) * 0.5 * eta2gaga();
968 
969  // Widths of decays Higgs -> Z0 + gamma0.
970  else if (id1Abs == 23 && id2Abs == 22)
971  widNow = preFac * pow2(alpEM / M_PI) * pow3(ps) * eta2gaZ();
972 
973  // Widths of decays Higgs (h0, H0) -> Z0 + Z0.
974  else if (id1Abs == 23 && id2Abs == 23) {
975  // If Higgs heavy use on-shell expression, else interpolation in table
976  if (mHat > 3. * mZ) kinFac = (1. - 4. * mr1 + 12. * mr1 * mr1) * ps;
977  else if (mHat > mLowZ) {
978  double xTab = (mHat - mLowZ) / mStepZ;
979  int iTab = max( 0, min( 99, int(xTab) ) );
980  kinFac = kinFacZ[iTab]
981  * pow( kinFacZ[iTab + 1] / kinFacZ[iTab], xTab - iTab );
982  }
983  else kinFac = 0.;
984  // Prefactor, normally rescaled to mRes^2 * mHat rather than mHat^3.
985  widNow = 0.25 * preFac * pow2(coup2Z) * kinFac;
986  if (!useCubicWidth) widNow *= pow2(mRes / mHat);
987  }
988 
989  // Widths of decays Higgs (h0, H0) -> W+ + W-.
990  else if (id1Abs == 24 && id2Abs == 24) {
991  // If Higgs heavy use on-shell expression, else interpolation in table.
992  if (mHat > 3. * mW) kinFac = (1. - 4. * mr1 + 12. * mr1 * mr1) * ps;
993  else if (mHat > mLowW) {
994  double xTab = (mHat - mLowW) / mStepW;
995  int iTab = max( 0, min( 99, int(xTab) ) );
996  kinFac = kinFacW[iTab]
997  * pow( kinFacW[iTab + 1] / kinFacW[iTab], xTab - iTab);
998  }
999  else kinFac = 0.;
1000  // Prefactor, normally rescaled to mRes^2 * mHat rather than mHat^3.
1001  widNow = 0.5 * preFac * pow2(coup2W) * kinFac;
1002  if (!useCubicWidth) widNow *= pow2(mRes / mHat);
1003  }
1004 
1005  // Widths of decays Higgs (H0) -> h0 + h0.
1006  else if (id1Abs == 25 && id2Abs == 25)
1007  widNow = 0.25 * preFac * pow4(mZ / mHat) * ps * pow2(coup2H1H1);
1008 
1009  // Widths of decays Higgs (H0) -> A0 + A0.
1010  else if (id1Abs == 36 && id2Abs == 36)
1011  widNow = 0.5 * preFac * pow4(mZ / mHat) * ps * pow2(coup2A3A3);
1012 
1013  // Widths of decays Higgs (A0) -> h0 + Z0.
1014  else if (id1Abs == 25 && id2Abs == 23)
1015  widNow = 0.5 * preFac * pow3(ps) * pow2(coup2H1Z);
1016 
1017  // Widths of decays Higgs (H0) -> A0 + Z0.
1018  else if (id1Abs == 36 && id2Abs == 23)
1019  widNow = 0.5 * preFac * pow3(ps) * pow2(coup2A3Z);
1020 
1021  // Widths of decays Higgs (H0) -> A0 + h0.
1022  else if (id1Abs == 36 && id2Abs == 25)
1023  widNow = 0.25 * preFac * pow4(mZ / mHat) * ps * pow2(coup2A3H1);
1024 
1025  // Widths of decays Higgs -> H+- + W-+.
1026  else if (id1Abs == 37 && id2Abs == 24)
1027  widNow = 0.5 * preFac * pow3(ps) * pow2(coup2HchgW);
1028 
1029 }
1030 
1031 //--------------------------------------------------------------------------
1032 
1033 // Sum up quark loop contributions in Higgs -> g + g.
1034 // Note: running quark masses are used, unlike Pythia6 (not negligible shift).
1035 
1036 double ResonanceH::eta2gg() {
1037 
1038  // Initial values.
1039  complex eta = complex(0., 0.);
1040  double mLoop, epsilon, root, rootLog;
1041  complex phi, etaNow;
1042 
1043  // Loop over s, c, b, t quark flavours.
1044  for (int idNow = 3; idNow < 7; ++idNow) {
1045  mLoop = (useRunLoopMass) ? particleDataPtr->mRun(idNow, mHat)
1046  : particleDataPtr->m0(idNow);
1047  epsilon = pow2(2. * mLoop / mHat);
1048 
1049  // Value of loop integral.
1050  if (epsilon <= 1.) {
1051  root = sqrt(1. - epsilon);
1052  rootLog = (epsilon < 1e-4) ? log(4. / epsilon - 2.)
1053  : log( (1. + root) / (1. - root) );
1054  phi = complex( -0.25 * (pow2(rootLog) - pow2(M_PI)),
1055  0.5 * M_PI * rootLog );
1056  }
1057  else phi = complex( pow2( asin(1. / sqrt(epsilon)) ), 0.);
1058 
1059  // Factors that depend on Higgs and flavour type.
1060  if (higgsType < 3) etaNow = -0.5 * epsilon
1061  * (complex(1., 0.) + (1. - epsilon) * phi);
1062  else etaNow = -0.5 * epsilon * phi;
1063  if (idNow%2 == 1) etaNow *= coup2d;
1064  else etaNow *= coup2u;
1065 
1066  // Sum up contribution and return square of absolute value.
1067  eta += etaNow;
1068  }
1069  return (pow2(eta.real()) + pow2(eta.imag()));
1070 
1071 }
1072 
1073 //--------------------------------------------------------------------------
1074 
1075 // Sum up quark, lepton, W+- and (for BSM) H+- loop contributions
1076 // in Higgs -> gamma + gamma.
1077 
1078 double ResonanceH::eta2gaga() {
1079 
1080  // Initial values.
1081  complex eta = complex(0., 0.);
1082  int idNow;
1083  double ef, mLoop, epsilon, root, rootLog;
1084  complex phi, etaNow;
1085 
1086  // Loop over s, c, b, t, mu, tau, W+-, H+- flavours.
1087  for (int idLoop = 0; idLoop < 8; ++idLoop) {
1088  if (idLoop < 4) idNow = idLoop + 3;
1089  else if (idLoop < 6) idNow = 2 * idLoop + 5;
1090  else if (idLoop < 7) idNow = 24;
1091  else idNow = 37;
1092  if (idNow == 37 && higgsType == 0) continue;
1093 
1094  // Charge and loop integral parameter.
1095  ef = (idNow < 20) ? couplingsPtr->ef(idNow) : 1.;
1096  mLoop = (useRunLoopMass) ? particleDataPtr->mRun(idNow, mHat)
1097  : particleDataPtr->m0(idNow);
1098  epsilon = pow2(2. * mLoop / mHat);
1099 
1100  // Value of loop integral.
1101  if (epsilon <= 1.) {
1102  root = sqrt(1. - epsilon);
1103  rootLog = (epsilon < 1e-4) ? log(4. / epsilon - 2.)
1104  : log( (1. + root) / (1. - root) );
1105  phi = complex( -0.25 * (pow2(rootLog) - pow2(M_PI)),
1106  0.5 * M_PI * rootLog );
1107  }
1108  else phi = complex( pow2( asin(1. / sqrt(epsilon)) ), 0.);
1109 
1110  // Expressions for quarks and leptons that depend on Higgs type.
1111  if (idNow < 17) {
1112  if (higgsType < 3) etaNow = -0.5 * epsilon
1113  * (complex(1., 0.) + (1. - epsilon) * phi);
1114  else etaNow = -0.5 * epsilon * phi;
1115  if (idNow < 7 && idNow%2 == 1) etaNow *= 3. * pow2(ef) * coup2d;
1116  else if (idNow < 7 ) etaNow *= 3. * pow2(ef) * coup2u;
1117  else etaNow *= pow2(ef) * coup2l;
1118  }
1119 
1120  // Expression for W+-.
1121  else if (idNow == 24) etaNow = (complex(0.5 + 0.75 * epsilon, 0.)
1122  + 0.75 * epsilon * (2. - epsilon) * phi) * coup2W;
1123 
1124  // Expression for H+-.
1125  else etaNow = (complex(epsilon, 0.) - epsilon * epsilon * phi)
1126  * pow2(mW / mHchg) * coup2Hchg;
1127 
1128  // Sum up contribution and return square of absolute value.
1129  eta += etaNow;
1130  }
1131  return (pow2(eta.real()) + pow2(eta.imag()));
1132 
1133 }
1134 
1135 //--------------------------------------------------------------------------
1136 
1137 // Sum up quark, lepton, W+- and (for BSM) H+- loop contributions
1138 // in Higgs -> gamma + Z0.
1139 
1140 double ResonanceH::eta2gaZ() {
1141 
1142  // Initial values.
1143  complex eta = complex(0., 0.);
1144  int idNow;
1145  double ef, vf, mLoop, epsilon, epsPrime, root, rootLog, asinEps;
1146  complex phi, psi, phiPrime, psiPrime, fXY, f1, etaNow;
1147 
1148  // Loop over s, c, b, t, mu , tau, W+-, H+- flavours.
1149  for (int idLoop = 0; idLoop < 7; ++idLoop) {
1150  if (idLoop < 4) idNow = idLoop + 3;
1151  else if (idLoop < 6) idNow = 2 * idLoop + 5;
1152  else if (idLoop < 7) idNow = 24;
1153  else idNow = 37;
1154 
1155  // Electroweak charges and loop integral parameters.
1156  ef = (idNow < 20) ? couplingsPtr->ef(idNow) : 1.;
1157  vf = (idNow < 20) ? couplingsPtr->vf(idNow) : 0.;
1158  mLoop = (useRunLoopMass) ? particleDataPtr->mRun(idNow, mHat)
1159  : particleDataPtr->m0(idNow);
1160  epsilon = pow2(2. * mLoop / mHat);
1161  epsPrime = pow2(2. * mLoop / mZ);
1162 
1163  // Value of loop integral for epsilon = 4 m^2 / sHat.
1164  if (epsilon <= 1.) {
1165  root = sqrt(1. - epsilon);
1166  rootLog = (epsilon < 1e-4) ? log(4. / epsilon - 2.)
1167  : log( (1. + root) / (1. - root) );
1168  phi = complex( -0.25 * (pow2(rootLog) - pow2(M_PI)),
1169  0.5 * M_PI * rootLog );
1170  psi = 0.5 * root * complex( rootLog, -M_PI);
1171  } else {
1172  asinEps = asin(1. / sqrt(epsilon));
1173  phi = complex( pow2(asinEps), 0.);
1174  psi = complex( sqrt(epsilon - 1.) * asinEps, 0.);
1175  }
1176 
1177  // Value of loop integral for epsilonPrime = 4 m^2 / m_Z^2.
1178  if (epsPrime <= 1.) {
1179  root = sqrt(1. - epsPrime);
1180  rootLog = (epsPrime < 1e-4) ? log(4. / epsPrime - 2.)
1181  : log( (1. + root) / (1. - root) );
1182  phiPrime = complex( -0.25 * (pow2(rootLog) - pow2(M_PI)),
1183  0.5 * M_PI * rootLog );
1184  psiPrime = 0.5 * root * complex( rootLog, -M_PI);
1185  } else {
1186  asinEps = asin(1. / sqrt(epsPrime));
1187  phiPrime = complex( pow2(asinEps), 0.);
1188  psiPrime = complex( sqrt(epsPrime - 1.) * asinEps, 0.);
1189  }
1190 
1191  // Combine the two loop integrals.
1192  fXY = (epsilon * epsPrime / (8. * pow2(epsilon - epsPrime)))
1193  * ( complex(epsilon - epsPrime, 0)
1194  + epsilon * epsPrime * (phi - phiPrime)
1195  + 2. * epsilon * (psi - psiPrime) );
1196  f1 = - (epsilon * epsPrime / (2. * (epsilon - epsPrime)))
1197  * (phi - phiPrime);
1198 
1199  // Expressions for quarks and leptons that depend on Higgs type.
1200  if (idNow < 17) {
1201  etaNow = (higgsType < 3) ? -fXY + 0.25 * f1 : 0.25 * f1;
1202  if (idNow < 7 && idNow%2 == 1) etaNow *= 3. * ef * vf * coup2d;
1203  else if (idNow < 7) etaNow *= 3. * ef * vf * coup2u;
1204  else etaNow *= ef * vf * coup2l;
1205 
1206  // Expression for W+-.
1207  } else if (idNow == 24) {
1208  double coef1 = 3. - sin2tW / cos2tW;
1209  double coefXY = (1. + 2. / epsilon) * sin2tW / cos2tW
1210  - (5. + 2. / epsilon);
1211  etaNow = -cos2tW * (coef1 * f1 + coefXY * fXY) * coup2W;
1212 
1213  // Expression for H+-.
1214  } else etaNow = (1. - 2. * sin2tW) * fXY * pow2(mW / mHchg)
1215  * coup2Hchg;
1216 
1217  // Sum up contribution and return square of absolute value.
1218  eta += etaNow;
1219  }
1220  return ( (pow2(eta.real()) + pow2(eta.imag())) / (sin2tW * cos2tW) );
1221 
1222 }
1223 
1224 //==========================================================================
1225 
1226 // The ResonanceHchg class.
1227 // Derived class for H+- properties.
1228 
1229 //--------------------------------------------------------------------------
1230 
1231 // Initialize constants.
1232 
1233 void ResonanceHchg::initConstants() {
1234 
1235  // Locally stored properties and couplings.
1236  useCubicWidth = settingsPtr->flag("Higgs:cubicWidth");
1237  thetaWRat = 1. / (8. * couplingsPtr->sin2thetaW());
1238  mW = particleDataPtr->m0(24);
1239  tanBeta = settingsPtr->parm("HiggsHchg:tanBeta");
1240  tan2Beta = tanBeta * tanBeta;
1241  coup2H1W = settingsPtr->parm("HiggsHchg:coup2H1W");
1242 
1243 }
1244 
1245 //--------------------------------------------------------------------------
1246 
1247 // Calculate various common prefactors for the current mass.
1248 
1249 void ResonanceHchg::calcPreFac(bool) {
1250 
1251  // Common coupling factors.
1252  alpEM = couplingsPtr->alphaEM(mHat * mHat);
1253  alpS = couplingsPtr->alphaS(mHat * mHat);
1254  colQ = 3. * (1. + alpS / M_PI);
1255  preFac = alpEM * thetaWRat * pow3(mHat) / pow2(mW);
1256 
1257 }
1258 
1259 //--------------------------------------------------------------------------
1260 
1261 // Calculate width for currently considered channel.
1262 
1263 void ResonanceHchg::calcWidth(bool) {
1264 
1265  // Check that above threshold.
1266  if (ps == 0.) return;
1267 
1268  // H+- decay to fermions involves running masses.
1269  if (id1Abs < 17 && (id1Abs < 7 || id1Abs > 10)) {
1270  double mRun1 = particleDataPtr->mRun(id1Abs, mHat);
1271  double mRun2 = particleDataPtr->mRun(id2Abs, mHat);
1272  double mrRunDn = pow2(mRun1 / mHat);
1273  double mrRunUp = pow2(mRun2 / mHat);
1274  if (id1Abs%2 == 0) swap( mrRunDn, mrRunUp);
1275 
1276  // Width to fermions: couplings, kinematics, colour factor.
1277  widNow = preFac * max( 0., (mrRunDn * tan2Beta + mrRunUp / tan2Beta)
1278  * (1. - mrRunDn - mrRunUp) - 4. *mrRunDn * mrRunUp ) * ps;
1279  if (id1Abs < 7) widNow *= colQ;
1280  }
1281 
1282  // H+- decay to h0 + W+-.
1283  else if (id1Abs == 25 && id2Abs == 24)
1284  widNow = 0.5 * preFac * pow3(ps) * pow2(coup2H1W);
1285 
1286 }
1287 
1288 //==========================================================================
1289 
1290 // The ResonanceZprime class.
1291 // Derived class for gamma*/Z0/Z'^0 properties.
1292 
1293 //--------------------------------------------------------------------------
1294 
1295 // Initialize constants.
1296 
1297 void ResonanceZprime::initConstants() {
1298 
1299  // Locally stored properties and couplings.
1300  gmZmode = settingsPtr->mode("Zprime:gmZmode");
1301  sin2tW = couplingsPtr->sin2thetaW();
1302  cos2tW = 1. - sin2tW;
1303  thetaWRat = 1. / (16. * sin2tW * cos2tW);
1304 
1305  // Properties of Z resonance.
1306  mZ = particleDataPtr->m0(23);
1307  GammaZ = particleDataPtr->mWidth(23);
1308  m2Z = mZ*mZ;
1309  GamMRatZ = GammaZ / mZ;
1310 
1311  // Ensure that arrays initially empty.
1312  for (int i = 0; i < 20; ++i) afZp[i] = 0.;
1313  for (int i = 0; i < 20; ++i) vfZp[i] = 0.;
1314 
1315  // Store first-generation axial and vector couplings.
1316  afZp[1] = settingsPtr->parm("Zprime:ad");
1317  afZp[2] = settingsPtr->parm("Zprime:au");
1318  afZp[11] = settingsPtr->parm("Zprime:ae");
1319  afZp[12] = settingsPtr->parm("Zprime:anue");
1320  vfZp[1] = settingsPtr->parm("Zprime:vd");
1321  vfZp[2] = settingsPtr->parm("Zprime:vu");
1322  vfZp[11] = settingsPtr->parm("Zprime:ve");
1323  vfZp[12] = settingsPtr->parm("Zprime:vnue");
1324 
1325  // Second and third generation could be carbon copy of this...
1326  if (settingsPtr->flag("Zprime:universality")) {
1327  for (int i = 3; i <= 6; ++i) {
1328  afZp[i] = afZp[i-2];
1329  vfZp[i] = vfZp[i-2];
1330  afZp[i+10] = afZp[i+8];
1331  vfZp[i+10] = vfZp[i+8];
1332  }
1333 
1334  // ... or could have different couplings.
1335  } else {
1336  afZp[3] = settingsPtr->parm("Zprime:as");
1337  afZp[4] = settingsPtr->parm("Zprime:ac");
1338  afZp[5] = settingsPtr->parm("Zprime:ab");
1339  afZp[6] = settingsPtr->parm("Zprime:at");
1340  afZp[13] = settingsPtr->parm("Zprime:amu");
1341  afZp[14] = settingsPtr->parm("Zprime:anumu");
1342  afZp[15] = settingsPtr->parm("Zprime:atau");
1343  afZp[16] = settingsPtr->parm("Zprime:anutau");
1344  vfZp[3] = settingsPtr->parm("Zprime:vs");
1345  vfZp[4] = settingsPtr->parm("Zprime:vc");
1346  vfZp[5] = settingsPtr->parm("Zprime:vb");
1347  vfZp[6] = settingsPtr->parm("Zprime:vt");
1348  vfZp[13] = settingsPtr->parm("Zprime:vmu");
1349  vfZp[14] = settingsPtr->parm("Zprime:vnumu");
1350  vfZp[15] = settingsPtr->parm("Zprime:vtau");
1351  vfZp[16] = settingsPtr->parm("Zprime:vnutau");
1352  }
1353 
1354  // Coupling for Z' -> W+ W-.
1355  coupZpWW = settingsPtr->parm("Zprime:coup2WW");
1356 
1357 }
1358 
1359 //--------------------------------------------------------------------------
1360 
1361 // Calculate various common prefactors for the current mass.
1362 
1363 void ResonanceZprime::calcPreFac(bool calledFromInit) {
1364 
1365  // Common coupling factors.
1366  alpEM = couplingsPtr->alphaEM(mHat * mHat);
1367  alpS = couplingsPtr->alphaS(mHat * mHat);
1368  colQ = 3. * (1. + alpS / M_PI);
1369  preFac = alpEM * thetaWRat * mHat / 3.;
1370 
1371  // When call for incoming flavour need to consider gamma*/Z0 mix.
1372  if (!calledFromInit) {
1373 
1374  // Couplings when an incoming fermion is specified; elso only pure Z'0.
1375  ei2 = 0.;
1376  eivi = 0.;
1377  vai2 = 0.;
1378  eivpi = 0.;
1379  vaivapi = 0.,
1380  vapi2 = 1.;
1381  int idInFlavAbs = abs(idInFlav);
1382  if (idInFlavAbs > 0 && idInFlavAbs < 19) {
1383  double ei = couplingsPtr->ef(idInFlavAbs);
1384  double ai = couplingsPtr->af(idInFlavAbs);
1385  double vi = couplingsPtr->vf(idInFlavAbs);
1386  double api = afZp[idInFlavAbs];
1387  double vpi = vfZp[idInFlavAbs];
1388  ei2 = ei * ei;
1389  eivi = ei * vi;
1390  vai2 = vi * vi + ai * ai;
1391  eivpi = ei * vpi;
1392  vaivapi = vi * vpi + ai * api;;
1393  vapi2 = vpi * vpi + api * api;
1394  }
1395 
1396  // Calculate prefactors for gamma/interference/Z0 terms.
1397  double sH = mHat * mHat;
1398  double propZ = sH / ( pow2(sH - m2Z) + pow2(sH * GamMRatZ) );
1399  double propZp = sH / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1400  gamNorm = ei2;
1401  gamZNorm = 2. * eivi * thetaWRat * (sH - m2Z) * propZ;
1402  ZNorm = vai2 * pow2(thetaWRat) * sH * propZ;
1403  gamZpNorm = 2. * eivpi * thetaWRat * (sH - m2Res) * propZp;
1404  ZZpNorm = 2. * vaivapi * pow2(thetaWRat) * ((sH - m2Res) * (sH - m2Z)
1405  + sH * GamMRat * sH * GamMRatZ) * propZ * propZp;
1406  ZpNorm = vapi2 * pow2(thetaWRat) * sH * propZp;
1407 
1408  // Rescale Z0 height normalization to compensate for a width one??
1409  //if (doForceWidth) {
1410  // intNorm *= forceFactor;
1411  // resNorm *= forceFactor;
1412  //}
1413 
1414  // Optionally only keep some of gamma*, Z0 and Z' terms.
1415  if (gmZmode == 1) {gamZNorm = 0; ZNorm = 0.; gamZpNorm = 0.;
1416  ZZpNorm = 0.; ZpNorm = 0.;}
1417  if (gmZmode == 2) {gamNorm = 0.; gamZNorm = 0.; gamZpNorm = 0.;
1418  ZZpNorm = 0.; ZpNorm = 0.;}
1419  if (gmZmode == 3) {gamNorm = 0.; gamZNorm = 0.; ZNorm = 0.;
1420  gamZpNorm = 0.; ZZpNorm = 0.;}
1421  if (gmZmode == 4) {gamZpNorm = 0.; ZZpNorm = 0.; ZpNorm = 0.;}
1422  if (gmZmode == 5) {gamZNorm = 0.; ZNorm = 0.; ZZpNorm = 0.;}
1423  if (gmZmode == 6) {gamNorm = 0.; gamZNorm = 0.; gamZpNorm = 0.;}
1424  }
1425 
1426 }
1427 
1428 //--------------------------------------------------------------------------
1429 
1430 // Calculate width for currently considered channel.
1431 
1432 void ResonanceZprime::calcWidth(bool calledFromInit) {
1433 
1434  // Check that above threshold.
1435  if (ps == 0.) return;
1436 
1437  // At initialization only the pure Z'0 should be considered.
1438  if (calledFromInit) {
1439 
1440  // Contributions from three fermion generations.
1441  if ( id1Abs < 7 || (id1Abs > 10 && id1Abs < 17) ) {
1442  double apf = afZp[id1Abs];
1443  double vpf = vfZp[id1Abs];
1444  widNow = preFac * ps * (vpf*vpf * (1. + 2. * mr1)
1445  + apf*apf * ps*ps);
1446  if (id1Abs < 7) widNow *= colQ;
1447 
1448  // Contribution from Z'0 -> W^+ W^-.
1449  } else if (id1Abs == 24) {
1450  widNow = preFac * pow2(coupZpWW * cos2tW) * pow3(ps)
1451  * (1. + mr1*mr1 + mr2*mr2 + 10. * (mr1 + mr2 + mr1 * mr2));
1452  }
1453  }
1454 
1455  // When call for incoming flavour need to consider full mix.
1456  else {
1457 
1458  // Contributions from three fermion generations.
1459  if ( id1Abs < 7 || (id1Abs > 10 && id1Abs < 17) ) {
1460 
1461  // Couplings of gamma^*/Z^0/Z'^0 to final flavour
1462  double ef = couplingsPtr->ef(id1Abs);
1463  double af = couplingsPtr->af(id1Abs);
1464  double vf = couplingsPtr->vf(id1Abs);
1465  double apf = afZp[id1Abs];
1466  double vpf = vfZp[id1Abs];
1467 
1468  // Combine couplings with kinematical factors.
1469  double kinFacA = pow3(ps);
1470  double kinFacV = ps * (1. + 2. * mr1);
1471  double ef2 = ef * ef * kinFacV;
1472  double efvf = ef * vf * kinFacV;
1473  double vaf2 = vf * vf * kinFacV + af * af * kinFacA;
1474  double efvpf = ef * vpf * kinFacV;
1475  double vafvapf = vf * vpf * kinFacV + af * apf * kinFacA;
1476  double vapf2 = vpf * vpf * kinFacV + apf * apf * kinFacA;
1477 
1478  // Relative outwidths: combine instate, propagator and outstate.
1479  widNow = gamNorm * ef2 + gamZNorm * efvf + ZNorm * vaf2
1480  + gamZpNorm * efvpf + ZZpNorm * vafvapf + ZpNorm * vapf2;
1481  if (id1Abs < 7) widNow *= colQ;
1482 
1483  // Contribution from Z'0 -> W^+ W^-.
1484  } else if (id1Abs == 24) {
1485  widNow = ZpNorm * pow2(coupZpWW * cos2tW) * pow3(ps)
1486  * (1. + mr1*mr1 + mr2*mr2 + 10. * (mr1 + mr2 + mr1 * mr2));
1487  }
1488  }
1489 
1490 }
1491 
1492 //==========================================================================
1493 
1494 // The ResonanceWprime class.
1495 // Derived class for W'+- properties.
1496 
1497 //--------------------------------------------------------------------------
1498 
1499 // Initialize constants.
1500 
1501 void ResonanceWprime::initConstants() {
1502 
1503  // Locally stored properties and couplings.
1504  thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());
1505  cos2tW = couplingsPtr->cos2thetaW();
1506 
1507  // Axial and vector couplings of fermions.
1508  aqWp = settingsPtr->parm("Wprime:aq");
1509  vqWp = settingsPtr->parm("Wprime:vq");
1510  alWp = settingsPtr->parm("Wprime:al");
1511  vlWp = settingsPtr->parm("Wprime:vl");
1512 
1513  // Coupling for W' -> W Z.
1514  coupWpWZ = settingsPtr->parm("Wprime:coup2WZ");
1515 
1516 }
1517 
1518 //--------------------------------------------------------------------------
1519 
1520 // Calculate various common prefactors for the current mass.
1521 
1522 void ResonanceWprime::calcPreFac(bool) {
1523 
1524  // Common coupling factors.
1525  alpEM = couplingsPtr->alphaEM(mHat * mHat);
1526  alpS = couplingsPtr->alphaS(mHat * mHat);
1527  colQ = 3. * (1. + alpS / M_PI);
1528  preFac = alpEM * thetaWRat * mHat;
1529 
1530 }
1531 
1532 //--------------------------------------------------------------------------
1533 
1534 // Calculate width for currently considered channel.
1535 
1536 void ResonanceWprime::calcWidth(bool) {
1537 
1538  // Check that above threshold.
1539  if (ps == 0.) return;
1540 
1541  // Decay to quarks involves colour factor and CKM matrix.
1542  if (id1Abs > 0 && id1Abs < 7) widNow
1543  = preFac * ps * 0.5 * ((vqWp*vqWp + aqWp * aqWp)
1544  + 6. * (vqWp*vqWp - aqWp * aqWp) * sqrt(mr1 *mr2))
1545  * (1. - 0.5 * (mr1 + mr2) - 0.5 * pow2(mr1 - mr2))
1546  * colQ * couplingsPtr->V2CKMid(id1Abs, id2Abs);
1547 
1548  // Decay to leptons simpler.
1549  else if (id1Abs > 10 && id1Abs < 17) widNow
1550  = preFac * ps * 0.5 * ((vlWp*vqWp + alWp * aqWp)
1551  + 6. * (vlWp*vqWp - alWp * aqWp) * sqrt(mr1 *mr2))
1552  * (1. - 0.5 * (mr1 + mr2) - 0.5 * pow2(mr1 - mr2));
1553 
1554  // Decay to W^+- Z^0.
1555  else if (id1Abs == 24 && id2Abs == 23) widNow
1556  = preFac * 0.25 * pow2(coupWpWZ) * cos2tW * (mr1 / mr2) * pow3(ps)
1557  * (1. + mr1*mr1 + mr2*mr2 + 10. * (mr1 + mr2 + mr1 * mr2));
1558 
1559 }
1560 
1561 //==========================================================================
1562 
1563 // The ResonanceRhorizontal class.
1564 // Derived class for R^0 (horizontal gauge boson) properties.
1565 
1566 //--------------------------------------------------------------------------
1567 
1568 // Initialize constants.
1569 
1570 void ResonanceRhorizontal::initConstants() {
1571 
1572  // Locally stored properties and couplings.
1573  thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());
1574 
1575 }
1576 
1577 //--------------------------------------------------------------------------
1578 
1579 // Calculate various common prefactors for the current mass.
1580 
1581 void ResonanceRhorizontal::calcPreFac(bool) {
1582 
1583  // Common coupling factors.
1584  alpEM = couplingsPtr->alphaEM(mHat * mHat);
1585  alpS = couplingsPtr->alphaS(mHat * mHat);
1586  colQ = 3. * (1. + alpS / M_PI);
1587  preFac = alpEM * thetaWRat * mHat;
1588 
1589 }
1590 
1591 //--------------------------------------------------------------------------
1592 
1593 // Calculate width for currently considered channel.
1594 
1595 void ResonanceRhorizontal::calcWidth(bool) {
1596 
1597  // Check that above threshold.
1598  if (ps == 0.) return;
1599 
1600  // R -> f fbar. Colour factor for quarks.
1601  widNow = preFac * ps * (2. - mr1 - mr2 - pow2(mr1 - mr2));
1602  if (id1Abs < 9) widNow *= colQ;
1603 
1604 }
1605 
1606 //==========================================================================
1607 
1608 // The ResonanceExcited class.
1609 // Derived class for excited-fermion properties.
1610 
1611 //--------------------------------------------------------------------------
1612 
1613 // Initialize constants.
1614 
1615 void ResonanceExcited::initConstants() {
1616 
1617  // Locally stored properties and couplings.
1618  Lambda = settingsPtr->parm("ExcitedFermion:Lambda");
1619  coupF = settingsPtr->parm("ExcitedFermion:coupF");
1620  coupFprime = settingsPtr->parm("ExcitedFermion:coupFprime");
1621  coupFcol = settingsPtr->parm("ExcitedFermion:coupFcol");
1622  sin2tW = couplingsPtr->sin2thetaW();
1623  cos2tW = 1. - sin2tW;
1624 
1625 }
1626 
1627 //--------------------------------------------------------------------------
1628 
1629 // Calculate various common prefactors for the current mass.
1630 
1631 void ResonanceExcited::calcPreFac(bool) {
1632 
1633  // Common coupling factors.
1634  alpEM = couplingsPtr->alphaEM(mHat * mHat);
1635  alpS = couplingsPtr->alphaS(mHat * mHat);
1636  preFac = pow3(mHat) / pow2(Lambda);
1637 
1638 }
1639 
1640 //--------------------------------------------------------------------------
1641 
1642 // Calculate width for currently considered channel.
1643 
1644 void ResonanceExcited::calcWidth(bool) {
1645 
1646  // Check that above threshold.
1647  if (ps == 0.) return;
1648 
1649  // f^* -> f g.
1650  if (id1Abs == 21) widNow = preFac * alpS * pow2(coupFcol) / 3.;
1651 
1652  // f^* -> f gamma.
1653  else if (id1Abs == 22) {
1654  double chgI3 = (id2Abs%2 == 0) ? 0.5 : -0.5;
1655  double chgY = (id2Abs < 9) ? 1. / 6. : -0.5;
1656  double chg = chgI3 * coupF + chgY * coupFprime;
1657  widNow = preFac * alpEM * pow2(chg) / 4.;
1658  }
1659 
1660  // f^* -> f Z^0.
1661  else if (id1Abs == 23) {
1662  double chgI3 = (id2Abs%2 == 0) ? 0.5 : -0.5;
1663  double chgY = (id2Abs < 9) ? 1. / 6. : -0.5;
1664  double chg = chgI3 * cos2tW * coupF - chgY * sin2tW * coupFprime;
1665  widNow = preFac * (alpEM * pow2(chg) / (8. * sin2tW * cos2tW))
1666  * ps*ps * (2. + mr1);
1667  }
1668 
1669  // f^* -> f' W^+-.
1670  else if (id1Abs == 24) widNow = preFac * (alpEM * pow2(coupF)
1671  / (16. * sin2tW)) * ps*ps * (2. + mr1);
1672 
1673 }
1674 
1675 //==========================================================================
1676 
1677 // The ResonanceGraviton class.
1678 // Derived class for excited Graviton properties.
1679 
1680 //--------------------------------------------------------------------------
1681 
1682 // Initialize constants.
1683 
1684 void ResonanceGraviton::initConstants() {
1685 
1686  // SMinBulk = off/on, use universal coupling (kappaMG)
1687  // or individual (Gxx) between graviton and SM particles.
1688  eDsmbulk = settingsPtr->flag("ExtraDimensionsG*:SMinBulk");
1689  eDvlvl = false;
1690  if (eDsmbulk) eDvlvl = settingsPtr->flag("ExtraDimensionsG*:VLVL");
1691  kappaMG = settingsPtr->parm("ExtraDimensionsG*:kappaMG");
1692  for (int i = 0; i < 27; ++i) eDcoupling[i] = 0.;
1693  double tmp_coup = settingsPtr->parm("ExtraDimensionsG*:Gqq");
1694  for (int i = 1; i <= 4; ++i) eDcoupling[i] = tmp_coup;
1695  eDcoupling[5] = settingsPtr->parm("ExtraDimensionsG*:Gbb");
1696  eDcoupling[6] = settingsPtr->parm("ExtraDimensionsG*:Gtt");
1697  tmp_coup = settingsPtr->parm("ExtraDimensionsG*:Gll");
1698  for (int i = 11; i <= 16; ++i) eDcoupling[i] = tmp_coup;
1699  eDcoupling[21] = settingsPtr->parm("ExtraDimensionsG*:Ggg");
1700  eDcoupling[22] = settingsPtr->parm("ExtraDimensionsG*:Ggmgm");
1701  eDcoupling[23] = settingsPtr->parm("ExtraDimensionsG*:GZZ");
1702  eDcoupling[24] = settingsPtr->parm("ExtraDimensionsG*:GWW");
1703  eDcoupling[25] = settingsPtr->parm("ExtraDimensionsG*:Ghh");
1704 
1705 }
1706 
1707 //--------------------------------------------------------------------------
1708 
1709 // Calculate various common prefactors for the current mass.
1710 
1711 void ResonanceGraviton::calcPreFac(bool) {
1712 
1713  // Common coupling factors.
1714  alpS = couplingsPtr->alphaS(mHat * mHat);
1715  colQ = 3. * (1. + alpS / M_PI);
1716  preFac = mHat / M_PI;
1717 
1718 }
1719 
1720 //--------------------------------------------------------------------------
1721 
1722 // Calculate width for currently considered channel.
1723 
1724 void ResonanceGraviton::calcWidth(bool) {
1725 
1726  // Check that above threshold.
1727  if (ps == 0.) return;
1728 
1729  // Widths to fermion pairs.
1730  if (id1Abs < 19) {
1731  widNow = preFac * pow3(ps) * (1. + 8. * mr1 / 3.) / 320.;
1732  if (id1Abs < 9) widNow *= colQ;
1733 
1734  // Widths to gluon and photon pair.
1735  } else if (id1Abs == 21) {
1736  widNow = preFac / 20.;
1737  } else if (id1Abs == 22) {
1738  widNow = preFac / 160.;
1739 
1740  // Widths to Z0 Z0 and W+ W- pair.
1741  } else if (id1Abs == 23 || id1Abs == 24) {
1742  // Longitudinal W/Z only.
1743  if (eDvlvl) {
1744  widNow = preFac * pow(ps,5) / 480.;
1745  // Transverse W/Z contributions as well.
1746  } else {
1747  widNow = preFac * ps * (13. / 12. + 14. * mr1 / 3. + 4. * mr1 * mr1)
1748  / 80.;
1749  }
1750  if (id1Abs == 23) widNow *= 0.5;
1751 
1752  // Widths to h h pair.
1753  } else if (id1Abs == 25) {
1754  widNow = preFac * pow(ps,5) / 960.;
1755  }
1756 
1757  // RS graviton coupling
1758  if (eDsmbulk) widNow *= 2. * pow2(eDcoupling[min( id1Abs, 26)] * mHat);
1759  else widNow *= pow2(kappaMG);
1760 
1761 }
1762 
1763 //==========================================================================
1764 
1765 // The ResonanceKKgluon class.
1766 // Derived class for excited kk-gluon properties.
1767 
1768 //--------------------------------------------------------------------------
1769 
1770 // Initialize constants.
1771 
1772 void ResonanceKKgluon::initConstants() {
1773 
1774  // KK-gluon gv/ga couplings and interference.
1775  for (int i = 0; i < 10; ++i) { eDgv[i] = 0.; eDga[i] = 0.; }
1776  double tmp_gL = settingsPtr->parm("ExtraDimensionsG*:KKgqL");
1777  double tmp_gR = settingsPtr->parm("ExtraDimensionsG*:KKgqR");
1778  for (int i = 1; i <= 4; ++i) {
1779  eDgv[i] = 0.5 * (tmp_gL + tmp_gR);
1780  eDga[i] = 0.5 * (tmp_gL - tmp_gR);
1781  }
1782  tmp_gL = settingsPtr->parm("ExtraDimensionsG*:KKgbL");
1783  tmp_gR = settingsPtr->parm("ExtraDimensionsG*:KKgbR");
1784  eDgv[5] = 0.5 * (tmp_gL + tmp_gR); eDga[5] = 0.5 * (tmp_gL - tmp_gR);
1785  tmp_gL = settingsPtr->parm("ExtraDimensionsG*:KKgtL");
1786  tmp_gR = settingsPtr->parm("ExtraDimensionsG*:KKgtR");
1787  eDgv[6] = 0.5 * (tmp_gL + tmp_gR); eDga[6] = 0.5 * (tmp_gL - tmp_gR);
1788  interfMode = settingsPtr->mode("ExtraDimensionsG*:KKintMode");
1789 
1790 }
1791 
1792 //--------------------------------------------------------------------------
1793 
1794 // Calculate various common prefactors for the current mass.
1795 
1796 void ResonanceKKgluon::calcPreFac(bool calledFromInit) {
1797 
1798  // Common coupling factors.
1799  alpS = couplingsPtr->alphaS(mHat * mHat);
1800  preFac = alpS * mHat / 6;
1801 
1802  // When call for incoming flavour need to consider g*/gKK mix.
1803  if (!calledFromInit) {
1804  // Calculate prefactors for g/interference/gKK terms.
1805  int idInFlavAbs = abs(idInFlav);
1806  double sH = mHat * mHat;
1807  normSM = 1;
1808  normInt = 2. * eDgv[min(idInFlavAbs, 9)] * sH * (sH - m2Res)
1809  / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1810  normKK = ( pow2(eDgv[min(idInFlavAbs, 9)])
1811  + pow2(eDga[min(idInFlavAbs, 9)]) ) * sH * sH
1812  / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1813 
1814  // Optionally only keep g* or gKK term.
1815  if (interfMode == 1) {normInt = 0.; normKK = 0.;}
1816  if (interfMode == 2) {normSM = 0.; normInt = 0.; normKK = 1.;}
1817  }
1818 
1819 }
1820 
1821 //--------------------------------------------------------------------------
1822 
1823 // Calculate width for currently considered channel.
1824 
1825 void ResonanceKKgluon::calcWidth(bool calledFromInit) {
1826 
1827  // Check that above threshold.
1828  if (ps == 0.) return;
1829 
1830  // Widths to quark pairs.
1831  if (id1Abs > 9) return;
1832 
1833  if (calledFromInit) {
1834  widNow = preFac * ps * (pow2(eDgv[min(id1Abs, 9)]) * (1. + 2.*mr1)
1835  + pow2(eDga[min(id1Abs, 9)]) * (1. - 4.*mr1) );
1836  } else {
1837  // Relative outwidths: combine instate, propagator and outstate.
1838  widNow = normSM * ps * (1. + 2. * mr1)
1839  + normInt * ps * eDgv[min(id1Abs, 9)] * (1. + 2. * mr1)
1840  + normKK * ps * (pow2(eDgv[min(id1Abs, 9)]) * (1. + 2.*mr1)
1841  + pow2(eDga[min(id1Abs, 9)]) * (1. - 4.*mr1) );
1842  widNow *= preFac;
1843  }
1844 
1845 }
1846 
1847 //==========================================================================
1848 
1849 // The ResonanceLeptoquark class.
1850 // Derived class for leptoquark properties.
1851 
1852 //--------------------------------------------------------------------------
1853 
1854 // Initialize constants.
1855 
1856 void ResonanceLeptoquark::initConstants() {
1857 
1858  // Locally stored properties and couplings.
1859  kCoup = settingsPtr->parm("LeptoQuark:kCoup");
1860 
1861  // Check that flavour info in decay channel is correctly set.
1862  int id1Now = particlePtr->channel(0).product(0);
1863  int id2Now = particlePtr->channel(0).product(1);
1864  if (id1Now < 1 || id1Now > 5) {
1865  infoPtr->errorMsg("Error in ResonanceLeptoquark::init:"
1866  " unallowed input quark flavour reset to u");
1867  id1Now = 2;
1868  particlePtr->channel(0).product(0, id1Now);
1869  }
1870  if (abs(id2Now) < 11 || abs(id2Now) > 16) {
1871  infoPtr->errorMsg("Error in ResonanceLeptoquark::init:"
1872  " unallowed input lepton flavour reset to e-");
1873  id2Now = 11;
1874  particlePtr->channel(0).product(1, id2Now);
1875  }
1876 
1877  // Set/overwrite charge and name of particle.
1878  bool changed = particlePtr->hasChanged();
1879  int chargeLQ = particleDataPtr->chargeType(id1Now)
1880  + particleDataPtr->chargeType(id2Now);
1881  particlePtr->setChargeType(chargeLQ);
1882  string nameLQ = "LQ_" + particleDataPtr->name(id1Now) + ","
1883  + particleDataPtr->name(id2Now);
1884  particlePtr->setNames(nameLQ, nameLQ + "bar");
1885  if (!changed) particlePtr->setHasChanged(false);
1886 
1887 }
1888 
1889 //--------------------------------------------------------------------------
1890 
1891 // Calculate various common prefactors for the current mass.
1892 
1893 void ResonanceLeptoquark::calcPreFac(bool) {
1894 
1895  // Common coupling factors.
1896  alpEM = couplingsPtr->alphaEM(mHat * mHat);
1897  preFac = 0.25 * alpEM * kCoup * mHat;
1898 
1899 }
1900 
1901 //--------------------------------------------------------------------------
1902 
1903 // Calculate width for currently considered channel.
1904 
1905 void ResonanceLeptoquark::calcWidth(bool) {
1906 
1907  // Check that above threshold.
1908  if (ps == 0.) return;
1909 
1910  // Width into lepton plus quark.
1911  if (id1Abs > 10 && id1Abs < 17 && id2Abs < 7) widNow = preFac * pow3(ps);
1912 
1913 }
1914 
1915 //==========================================================================
1916 
1917 // The ResonanceNuRight class.
1918 // Derived class for righthanded Majorana neutrino properties.
1919 
1920 //--------------------------------------------------------------------------
1921 
1922 // Initialize constants.
1923 
1924 void ResonanceNuRight::initConstants() {
1925 
1926  // Locally stored properties and couplings: righthanded W mass.
1927  thetaWRat = 1. / (768. * M_PI * pow2(couplingsPtr->sin2thetaW()));
1928  mWR = particleDataPtr->m0(9900024);
1929 
1930 }
1931 
1932 //--------------------------------------------------------------------------
1933 
1934 // Calculate various common prefactors for the current mass.
1935 
1936 void ResonanceNuRight::calcPreFac(bool) {
1937 
1938  // Common coupling factors.
1939  alpEM = couplingsPtr->alphaEM(mHat * mHat);
1940  alpS = couplingsPtr->alphaS(mHat * mHat);
1941  colQ = 3. * (1. + alpS / M_PI);
1942  preFac = pow2(alpEM) * thetaWRat * pow5(mHat) / pow4(max(mHat, mWR));
1943 
1944 }
1945 
1946 //--------------------------------------------------------------------------
1947 
1948 // Calculate width for currently considered channel.
1949 
1950 void ResonanceNuRight::calcWidth(bool) {
1951 
1952  // Check that above threshold.
1953  if (mHat < mf1 + mf2 + mf3 + MASSMARGIN) return;
1954 
1955  // Coupling part of widths to l- q qbar', l- l'+ nu_lR' and c.c.
1956  widNow = (id2Abs < 9 && id3Abs < 9)
1957  ? preFac * colQ * couplingsPtr->V2CKMid(id2, id3) : preFac;
1958 
1959  // Phase space corrections in decay. Must have y < 1.
1960  double x = (mf1 + mf2 + mf3) / mHat;
1961  double x2 = x * x;
1962  double fx = 1. - 8. * x2 + 8. * pow3(x2) - pow4(x2)
1963  - 24. * pow2(x2) * log(x);
1964  double y = min( 0.999, pow2(mHat / mWR) );
1965  double fy = ( 12. * (1. - y) * log(1. - y) + 12. * y - 6. * y*y
1966  - 2.* pow3(y) ) / pow4(y);
1967  widNow *= fx * fy;
1968 
1969 }
1970 
1971 //==========================================================================
1972 
1973 // The ResonanceZRight class.
1974 // Derived class for Z_R^0 properties.
1975 
1976 //--------------------------------------------------------------------------
1977 
1978 // Initialize constants.
1979 
1980 void ResonanceZRight::initConstants() {
1981 
1982  // Locally stored properties and couplings: righthanded W mass.
1983  sin2tW = couplingsPtr->sin2thetaW();
1984  thetaWRat = 1. / (48. * sin2tW * (1. - sin2tW) * (1. - 2. * sin2tW) );
1985 
1986 }
1987 
1988 //--------------------------------------------------------------------------
1989 
1990 // Calculate various common prefactors for the current mass.
1991 
1992 void ResonanceZRight::calcPreFac(bool) {
1993 
1994  // Common coupling factors.
1995  alpEM = couplingsPtr->alphaEM(mHat * mHat);
1996  alpS = couplingsPtr->alphaS(mHat * mHat);
1997  colQ = 3. * (1. + alpS / M_PI);
1998  preFac = alpEM * thetaWRat * mHat;
1999 
2000 }
2001 
2002 //--------------------------------------------------------------------------
2003 
2004 // Calculate width for currently considered channel.
2005 
2006 void ResonanceZRight::calcWidth(bool) {
2007 
2008  // Check that above threshold.
2009  if (ps == 0.) return;
2010 
2011  // Couplings to q qbar and l+ l-.
2012  double vf = 0.;
2013  double af = 0.;
2014  double symMaj = 1.;
2015  if (id1Abs < 9 && id1Abs%2 == 1) {
2016  af = -1. + 2. * sin2tW;
2017  vf = -1. + 4. * sin2tW / 3.;
2018  } else if (id1Abs < 9) {
2019  af = 1. - 2. * sin2tW;
2020  vf = 1. - 8. * sin2tW / 3.;
2021  } else if (id1Abs < 19 && id1Abs%2 == 1) {
2022  af = -1. + 2. * sin2tW;
2023  vf = -1. + 4. * sin2tW;
2024 
2025  // Couplings to nu_L nu_Lbar and nu_R nu_Rbar, both assumed Majoranas.
2026  } else if (id1Abs < 19) {
2027  af = -2. * sin2tW;
2028  vf = 0.;
2029  symMaj = 0.5;
2030  } else {
2031  af = 2. * (1. - sin2tW);
2032  vf = 0.;
2033  symMaj = 0.5;
2034  }
2035 
2036  // Width expression, including phase space and colour factor.
2037  widNow = preFac * (vf*vf * (1. + 2. * mr1) + af*af * ps*ps) * ps
2038  * symMaj;
2039  if (id1Abs < 9) widNow *= colQ;
2040 
2041 }
2042 
2043 //==========================================================================
2044 
2045 // The ResonanceWRight class.
2046 // Derived class for W_R+- properties.
2047 
2048 //--------------------------------------------------------------------------
2049 
2050 // Initialize constants.
2051 
2052 void ResonanceWRight::initConstants() {
2053 
2054  // Locally stored properties and couplings.
2055  thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());
2056 
2057 }
2058 
2059 //--------------------------------------------------------------------------
2060 
2061 // Calculate various common prefactors for the current mass.
2062 
2063 void ResonanceWRight::calcPreFac(bool) {
2064 
2065  // Common coupling factors.
2066  alpEM = couplingsPtr->alphaEM(mHat * mHat);
2067  alpS = couplingsPtr->alphaS(mHat * mHat);
2068  colQ = 3. * (1. + alpS / M_PI);
2069  preFac = alpEM * thetaWRat * mHat;
2070 
2071 }
2072 
2073 //--------------------------------------------------------------------------
2074 
2075 // Calculate width for currently considered channel.
2076 
2077 void ResonanceWRight::calcWidth(bool) {
2078 
2079  // Check that above threshold.
2080  if (ps == 0.) return;
2081 
2082  // Combine kinematics with colour factor and CKM couplings.
2083  widNow = preFac * (1. - 0.5 * (mr1 + mr2) - 0.5 * pow2(mr1 - mr2))
2084  * ps;
2085  if (id1Abs < 9) widNow *= colQ * couplingsPtr->V2CKMid(id1Abs, id2Abs);
2086 
2087 }
2088 
2089 //==========================================================================
2090 
2091 // The ResonanceHchgchgLeft class.
2092 // Derived class for H++/H-- (left) properties.
2093 
2094 //--------------------------------------------------------------------------
2095 
2096 // Initialize constants.
2097 
2098 void ResonanceHchgchgLeft::initConstants() {
2099 
2100  // Read in Yukawa matrix for couplings to a lepton pair.
2101  yukawa[1][1] = settingsPtr->parm("LeftRightSymmmetry:coupHee");
2102  yukawa[2][1] = settingsPtr->parm("LeftRightSymmmetry:coupHmue");
2103  yukawa[2][2] = settingsPtr->parm("LeftRightSymmmetry:coupHmumu");
2104  yukawa[3][1] = settingsPtr->parm("LeftRightSymmmetry:coupHtaue");
2105  yukawa[3][2] = settingsPtr->parm("LeftRightSymmmetry:coupHtaumu");
2106  yukawa[3][3] = settingsPtr->parm("LeftRightSymmmetry:coupHtautau");
2107 
2108  // Locally stored properties and couplings.
2109  gL = settingsPtr->parm("LeftRightSymmmetry:gL");
2110  vL = settingsPtr->parm("LeftRightSymmmetry:vL");
2111  mW = particleDataPtr->m0(24);
2112 
2113 }
2114 
2115 //--------------------------------------------------------------------------
2116 
2117 // Calculate various common prefactors for the current mass.
2118 
2119 void ResonanceHchgchgLeft::calcPreFac(bool) {
2120 
2121  // Common coupling factors.
2122  preFac = mHat / (8. * M_PI);
2123 
2124 }
2125 
2126 //--------------------------------------------------------------------------
2127 
2128 // Calculate width for currently considered channel.
2129 
2130 void ResonanceHchgchgLeft::calcWidth(bool) {
2131 
2132  // Check that above threshold.
2133  if (ps == 0.) return;
2134 
2135  // H++-- width to a pair of leptons. Combinatorial factor of 2.
2136  if (id1Abs < 17 && id2Abs < 17) {
2137  widNow = preFac * pow2(yukawa[(id1Abs-9)/2][(id2Abs-9)/2]) * ps;
2138  if (id2Abs != id1Abs) widNow *= 2.;
2139  }
2140 
2141  // H++-- width to a pair of lefthanded W's.
2142  else if (id1Abs == 24 && id2Abs == 24)
2143  widNow = preFac * 0.5 * pow2(gL*gL * vL / mW)
2144  * (3. * mr1 + 0.25 / mr1 - 1.) * ps;
2145 
2146 }
2147 
2148 //==========================================================================
2149 
2150 // The ResonanceHchgchgRight class.
2151 // Derived class for H++/H-- (right) properties.
2152 
2153 //--------------------------------------------------------------------------
2154 
2155 // Initialize constants.
2156 
2157 void ResonanceHchgchgRight::initConstants() {
2158 
2159  // Read in Yukawa matrix for couplings to a lepton pair.
2160  yukawa[1][1] = settingsPtr->parm("LeftRightSymmmetry:coupHee");
2161  yukawa[2][1] = settingsPtr->parm("LeftRightSymmmetry:coupHmue");
2162  yukawa[2][2] = settingsPtr->parm("LeftRightSymmmetry:coupHmumu");
2163  yukawa[3][1] = settingsPtr->parm("LeftRightSymmmetry:coupHtaue");
2164  yukawa[3][2] = settingsPtr->parm("LeftRightSymmmetry:coupHtaumu");
2165  yukawa[3][3] = settingsPtr->parm("LeftRightSymmmetry:coupHtautau");
2166 
2167  // Locally stored properties and couplings.
2168  idWR = 9000024;
2169  gR = settingsPtr->parm("LeftRightSymmmetry:gR");
2170 
2171 }
2172 
2173 //--------------------------------------------------------------------------
2174 
2175 // Calculate various common prefactors for the current mass.
2176 
2177 void ResonanceHchgchgRight::calcPreFac(bool) {
2178 
2179  // Common coupling factors.
2180  preFac = mHat / (8. * M_PI);
2181 
2182 }
2183 
2184 //--------------------------------------------------------------------------
2185 
2186 // Calculate width for currently considered channel.
2187 
2188 void ResonanceHchgchgRight::calcWidth(bool) {
2189 
2190  // Check that above threshold.
2191  if (ps == 0.) return;
2192 
2193  // H++-- width to a pair of leptons. Combinatorial factor of 2.
2194  if (id1Abs < 17 && id2Abs < 17) {
2195  widNow = preFac * pow2(yukawa[(id1Abs-9)/2][(id2Abs-9)/2]) * ps;
2196  if (id2Abs != id1Abs) widNow *= 2.;
2197  }
2198 
2199  // H++-- width to a pair of lefthanded W's.
2200  else if (id1Abs == idWR && id2Abs == idWR)
2201  widNow = preFac * pow2(yukawa[(id1Abs-9)/2][(id2Abs-9)/2]) * ps;
2202 
2203 }
2204 
2205 //==========================================================================
2206 
2207 } // end namespace Pythia8
Generated by   doxygen 1.8.5