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