Update to pythi8.170
[u/mrichter/AliRoot.git] / PYTHIA8 / pythia8170 / src / ResonanceWidths.cxx
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   // Extra coupling factors for t -> H+ + b.
724   tanBeta   = settingsPtr->parm("HiggsHchg:tanBeta");
725   tan2Beta  = tanBeta * tanBeta;
726   mbRun     = particleDataPtr->mRun( 5, particleDataPtr->m0(6) ); 
727  
728 }
729
730 //--------------------------------------------------------------------------
731
732 // Calculate various common prefactors for the current mass.
733
734 void ResonanceTop::calcPreFac(bool) {
735
736   // Common coupling factors.
737   alpEM     = couplingsPtr->alphaEM(mHat * mHat);
738   alpS      = couplingsPtr->alphaS(mHat * mHat);
739   colQ      = 1. - 2.5 * alpS / M_PI;
740   preFac    = alpEM * thetaWRat * pow3(mHat) / m2W;
741
742 }
743
744 //--------------------------------------------------------------------------
745
746 // Calculate width for currently considered channel.
747
748 void ResonanceTop::calcWidth(bool) {
749
750
751   // Check that above threshold. 
752   if (ps == 0.) return;
753
754   // Contributions from W + quark.
755   if (id1Abs == 24 && id2Abs < 6) { 
756     widNow  = preFac * ps 
757             * ( pow2(1. - mr2) + (1. + mr2) * mr1 - 2. * mr1 * mr1 );
758
759     // Combine with colour factor and CKM couplings.
760     widNow *= colQ * couplingsPtr->V2CKMid(6, id2Abs);
761
762   // Contributions from H+ + quark (so far only b).
763   } else if (id1Abs == 37 && id2Abs == 5) {     
764     widNow  = preFac * ps * ( (1. + mr2 - mr1) 
765             * (pow2(mbRun / mHat) * tan2Beta + 1. / tan2Beta) 
766             + 4. * mbRun * mf2 / pow2(mHat) ); 
767   }
768
769 }
770  
771 //==========================================================================
772
773 // The ResonanceFour class.
774 // Derived class for fourth-generation properties.
775
776 //--------------------------------------------------------------------------
777
778 // Initialize constants.
779
780 void ResonanceFour::initConstants() {
781
782   // Locally stored properties and couplings.
783   thetaWRat = 1. / (16. * couplingsPtr->sin2thetaW());
784   m2W       = pow2(particleDataPtr->m0(24));
785
786 }
787
788 //--------------------------------------------------------------------------
789
790 // Calculate various common prefactors for the current mass.
791
792 void ResonanceFour::calcPreFac(bool) {
793
794   // Common coupling factors.
795   alpEM     = couplingsPtr->alphaEM(mHat * mHat);
796   alpS      = couplingsPtr->alphaS(mHat * mHat);
797   colQ      = (idRes < 9) ? 1. - 2.5 * alpS / M_PI : 1.;
798   preFac    = alpEM * thetaWRat * pow3(mHat) / m2W;
799
800 }
801
802 //--------------------------------------------------------------------------
803
804 // Calculate width for currently considered channel.
805
806 void ResonanceFour::calcWidth(bool) {
807
808   // Only contributions from W + fermion.
809   if (id1Abs != 24 || id2Abs > 18) return; 
810
811   // Check that above threshold. Kinematical factor.
812   if (ps == 0.) return;
813   widNow    = preFac * ps 
814             * ( pow2(1. - mr2) + (1. + mr2) * mr1 - 2. * mr1 * mr1 );
815
816   // Combine with colour factor and CKM couplings.
817   if (idRes < 9) widNow *= colQ * couplingsPtr->V2CKMid(idRes, id2Abs);
818
819 }
820  
821 //==========================================================================
822
823 // The ResonanceH class.
824 // Derived class for SM and BSM Higgs properties.
825
826 //--------------------------------------------------------------------------
827
828 // Constants: could be changed here if desired, but normally should not.
829 // These are of technical nature, as described for each.
830
831 // Minimal mass for W, Z, top in integration over respective Breit-Wigner.
832 // Top constrainted by t -> W b decay, which is not seen in simple top BW.
833 const double ResonanceH::MASSMINWZ = 10.;
834 const double ResonanceH::MASSMINT  = 100.;
835
836 // Number of widths above threshold where B-W integration not needed.
837 const double ResonanceH::GAMMAMARGIN = 10.;
838
839 //--------------------------------------------------------------------------
840
841 // Initialize constants.
842
843 void ResonanceH::initConstants() {
844
845   // Locally stored properties and couplings.
846   useCubicWidth  = settingsPtr->flag("Higgs:cubicWidth");
847   useRunLoopMass = settingsPtr->flag("Higgs:runningLoopMass");
848   sin2tW         = couplingsPtr->sin2thetaW();
849   cos2tW         = 1. - sin2tW;
850   mT             = particleDataPtr->m0(6);
851   mZ             = particleDataPtr->m0(23);
852   mW             = particleDataPtr->m0(24);
853   mHchg          = particleDataPtr->m0(37);
854   GammaT         = particleDataPtr->mWidth(6);
855   GammaZ         = particleDataPtr->mWidth(23);
856   GammaW         = particleDataPtr->mWidth(24);
857
858   // Couplings to fermions, Z and W, depending on Higgs type.
859   coup2d         = 1.;
860   coup2u         = 1.;
861   coup2l         = 1.;
862   coup2Z         = 1.;
863   coup2W         = 1.;
864   coup2Hchg      = 0.;
865   coup2H1H1      = 0.;
866   coup2A3A3      = 0.;
867   coup2H1Z       = 0.;
868   coup2A3Z       = 0.;
869   coup2A3H1      = 0.;
870   coup2HchgW     = 0.;
871   if (higgsType == 1) {
872     coup2d       = settingsPtr->parm("HiggsH1:coup2d");
873     coup2u       = settingsPtr->parm("HiggsH1:coup2u");
874     coup2l       = settingsPtr->parm("HiggsH1:coup2l");
875     coup2Z       = settingsPtr->parm("HiggsH1:coup2Z");
876     coup2W       = settingsPtr->parm("HiggsH1:coup2W");
877     coup2Hchg    = settingsPtr->parm("HiggsH1:coup2Hchg");
878   } else if (higgsType == 2) {
879     coup2d       = settingsPtr->parm("HiggsH2:coup2d");
880     coup2u       = settingsPtr->parm("HiggsH2:coup2u");
881     coup2l       = settingsPtr->parm("HiggsH2:coup2l");
882     coup2Z       = settingsPtr->parm("HiggsH2:coup2Z");
883     coup2W       = settingsPtr->parm("HiggsH2:coup2W");
884     coup2Hchg    = settingsPtr->parm("HiggsH2:coup2Hchg");
885     coup2H1H1    = settingsPtr->parm("HiggsH2:coup2H1H1");
886     coup2A3A3    = settingsPtr->parm("HiggsH2:coup2A3A3");
887     coup2H1Z     = settingsPtr->parm("HiggsH2:coup2H1Z");
888     coup2A3Z     = settingsPtr->parm("HiggsA3:coup2H2Z");
889     coup2A3H1    = settingsPtr->parm("HiggsH2:coup2A3H1");
890     coup2HchgW   = settingsPtr->parm("HiggsH2:coup2HchgW");
891   } else if (higgsType == 3) {
892     coup2d       = settingsPtr->parm("HiggsA3:coup2d");
893     coup2u       = settingsPtr->parm("HiggsA3:coup2u");
894     coup2l       = settingsPtr->parm("HiggsA3:coup2l");
895     coup2Z       = settingsPtr->parm("HiggsA3:coup2Z");
896     coup2W       = settingsPtr->parm("HiggsA3:coup2W");
897     coup2Hchg    = settingsPtr->parm("HiggsA3:coup2Hchg");
898     coup2H1H1    = settingsPtr->parm("HiggsA3:coup2H1H1");
899     coup2H1Z     = settingsPtr->parm("HiggsA3:coup2H1Z");
900     coup2HchgW   = settingsPtr->parm("HiggsA3:coup2Hchg");
901   }
902
903   // Initialization of threshold kinematical factor by stepwise
904   // numerical integration of H -> t tbar, Z0 Z0 and W+ W-. 
905   int psModeT  = (higgsType < 3) ? 3 : 1;       
906   int psModeWZ = (higgsType < 3) ? 5 : 6; 
907   mLowT        = max( 2.02 * MASSMINT, 0.5 * mT);
908   mStepT       = 0.01 * (3. * mT - mLowT);    
909   mLowZ        = max( 2.02 * MASSMINWZ, 0.5 * mZ);
910   mStepZ       = 0.01 * (3. * mZ - mLowZ);    
911   mLowW        = max( 2.02 * MASSMINWZ, 0.5 * mW);
912   mStepW       = 0.01 * (3. * mW - mLowW);    
913   for (int i = 0; i <= 100; ++i) { 
914     kinFacT[i] = numInt2BW( mLowT + i * mStepT, 
915                  mT, GammaT, MASSMINT,  mT, GammaT, MASSMINT,  psModeT);
916     kinFacZ[i] = numInt2BW( mLowZ + i * mStepZ,
917                  mZ, GammaZ, MASSMINWZ, mZ, GammaZ, MASSMINWZ, psModeWZ);
918     kinFacW[i] = numInt2BW( mLowW + i * mStepW,
919                  mW, GammaW, MASSMINWZ, mW, GammaW, MASSMINWZ, psModeWZ);
920   }
921
922 }
923
924 //--------------------------------------------------------------------------
925
926 // Calculate various common prefactors for the current mass.
927
928 void ResonanceH::calcPreFac(bool) {
929
930   // Common coupling factors.
931   alpEM     = couplingsPtr->alphaEM(mHat * mHat);
932   alpS      = couplingsPtr->alphaS(mHat * mHat);
933   colQ      = 3. * (1. + alpS / M_PI);
934   preFac    = (alpEM / (8. * sin2tW)) * pow3(mHat) / pow2(mW); 
935
936 }
937
938 //--------------------------------------------------------------------------
939
940 // Calculate width for currently considered channel.
941
942 void ResonanceH::calcWidth(bool) {
943
944   // Widths of decays Higgs -> f + fbar.
945   if ( id2Abs == id1Abs && ( (id1Abs > 0 && id1Abs < 7) 
946     || (id1Abs > 10 && id1Abs < 17) ) ) {
947     kinFac = 0.;
948
949     // Check that above threshold (well above for top). Kinematical factor.
950     if ( (id1Abs != 6 && mHat > 2. * mf1 + MASSMARGIN) 
951       || (id1Abs == 6 && mHat > 3. * mT ) ) {
952       // A0 behaves like beta, h0 and H0 like beta**3.
953       kinFac = (higgsType < 3) ? pow3(ps) : ps;
954     }
955
956     // Top near or below threshold: interpolate in table.
957     else if (id1Abs == 6 && mHat > mLowT) {
958       double xTab = (mHat - mLowT) / mStepT;
959       int    iTab = max( 0, min( 99, int(xTab) ) );
960       kinFac      = kinFacT[iTab] 
961                   * pow( kinFacT[iTab + 1] / kinFacT[iTab], xTab - iTab);
962     }
963
964     // Coupling from mass and from BSM deviation from SM.
965     double coupFac = pow2(particleDataPtr->mRun(id1Abs, mHat) / mHat);
966     if (id1Abs < 7 && id1Abs%2 == 1) coupFac *= coup2d * coup2d;
967     else if (id1Abs < 7)             coupFac *= coup2u * coup2u;   
968     else                             coupFac *= coup2l * coup2l;
969
970     // Combine couplings and phase space with colour factor.
971     widNow = preFac * coupFac * kinFac; 
972     if (id1Abs < 7) widNow *= colQ;
973   }
974
975   // Widths of decays Higgs -> g + g.
976   else if (id1Abs == 21 && id2Abs == 21)  
977     widNow = preFac * pow2(alpS / M_PI) * eta2gg(); 
978
979   // Widths of decays Higgs -> gamma + gamma.
980   else if (id1Abs == 22 && id2Abs == 22)  
981     widNow = preFac * pow2(alpEM / M_PI) * 0.5 * eta2gaga(); 
982  
983   // Widths of decays Higgs -> Z0 + gamma0.
984   else if (id1Abs == 23 && id2Abs == 22) 
985     widNow = preFac * pow2(alpEM / M_PI) * pow3(ps) * eta2gaZ(); 
986     
987   // Widths of decays Higgs (h0, H0) -> Z0 + Z0.
988   else if (id1Abs == 23 && id2Abs == 23) {
989     // If Higgs heavy use on-shell expression, else interpolation in table
990     if (mHat > 3. * mZ) kinFac = (1.  - 4. * mr1 + 12. * mr1 * mr1) * ps;
991     else if (mHat > mLowZ) {
992       double xTab = (mHat - mLowZ) / mStepZ;
993       int    iTab = max( 0, min( 99, int(xTab) ) );
994       kinFac      = kinFacZ[iTab] 
995                   * pow( kinFacZ[iTab + 1] / kinFacZ[iTab], xTab - iTab );
996     }
997     else kinFac   = 0.;
998     // Prefactor, normally rescaled to mRes^2 * mHat rather than mHat^3.
999     widNow        = 0.25 * preFac * pow2(coup2Z) * kinFac;
1000     if (!useCubicWidth) widNow *= pow2(mRes / mHat);   
1001   }
1002  
1003   // Widths of decays Higgs (h0, H0) -> W+ + W-.
1004   else if (id1Abs == 24 && id2Abs == 24) {
1005     // If Higgs heavy use on-shell expression, else interpolation in table.
1006     if (mHat > 3. * mW) kinFac = (1.  - 4. * mr1 + 12. * mr1 * mr1) * ps;
1007     else if (mHat > mLowW) {
1008       double xTab = (mHat - mLowW) / mStepW;
1009       int    iTab = max( 0, min( 99, int(xTab) ) );
1010       kinFac      = kinFacW[iTab] 
1011                   * pow( kinFacW[iTab + 1] / kinFacW[iTab], xTab - iTab);
1012     }
1013     else kinFac   = 0.;
1014     // Prefactor, normally rescaled to mRes^2 * mHat rather than mHat^3.
1015     widNow        = 0.5 * preFac * pow2(coup2W) * kinFac;
1016     if (!useCubicWidth) widNow *= pow2(mRes / mHat);   
1017   }
1018  
1019   // Widths of decays Higgs (H0) -> h0 + h0.
1020   else if (id1Abs == 25 && id2Abs == 25) 
1021     widNow = 0.25 * preFac * pow4(mZ / mHat) * ps * pow2(coup2H1H1);
1022      
1023   // Widths of decays Higgs (H0) -> A0 + A0.
1024   else if (id1Abs == 36 && id2Abs == 36) 
1025     widNow = 0.5 * preFac * pow4(mZ / mHat) * ps * pow2(coup2A3A3);
1026  
1027   // Widths of decays Higgs (A0) -> h0 + Z0.
1028   else if (id1Abs == 25 && id2Abs == 23) 
1029     widNow = 0.5 * preFac * pow3(ps) * pow2(coup2H1Z);
1030  
1031   // Widths of decays Higgs (H0) -> A0 + Z0.
1032   else if (id1Abs == 36 && id2Abs == 23) 
1033     widNow = 0.5 * preFac * pow3(ps) * pow2(coup2A3Z);
1034    
1035   // Widths of decays Higgs (H0) -> A0 + h0.
1036   else if (id1Abs == 36 && id2Abs == 25) 
1037     widNow = 0.25 * preFac * pow4(mZ / mHat) * ps * pow2(coup2A3H1);
1038  
1039   // Widths of decays Higgs -> H+- + W-+.
1040   else if (id1Abs == 37 && id2Abs == 24) 
1041     widNow = 0.5 * preFac * pow3(ps) * pow2(coup2HchgW);
1042
1043 }
1044
1045 //--------------------------------------------------------------------------
1046
1047 // Sum up quark loop contributions in Higgs -> g + g.
1048 // Note: running quark masses are used, unlike Pythia6 (not negligible shift). 
1049
1050 double ResonanceH::eta2gg() {
1051
1052   // Initial values.
1053   complex eta = complex(0., 0.);
1054   double  mLoop, epsilon, root, rootLog;
1055   complex phi, etaNow;
1056
1057   // Loop over s, c, b, t quark flavours.
1058   for (int idNow = 3; idNow < 7; ++idNow) {
1059     mLoop   = (useRunLoopMass) ? particleDataPtr->mRun(idNow, mHat)
1060                                : particleDataPtr->m0(idNow);
1061     epsilon = pow2(2. * mLoop / mHat);
1062
1063     // Value of loop integral.
1064     if (epsilon <= 1.) {
1065       root    = sqrt(1. - epsilon);
1066       rootLog = (epsilon < 1e-4) ? log(4. / epsilon - 2.)
1067                 : log( (1. + root) / (1. - root) );
1068       phi     = complex( -0.25 * (pow2(rootLog) - pow2(M_PI)), 
1069                 0.5 * M_PI * rootLog );
1070     } 
1071     else phi  = complex( pow2( asin(1. / sqrt(epsilon)) ), 0.);
1072   
1073     // Factors that depend on Higgs and flavour type.
1074     if (higgsType < 3) etaNow = -0.5 * epsilon 
1075       * (complex(1., 0.) + (1. - epsilon) * phi);
1076     else etaNow = -0.5 * epsilon * phi;    
1077     if (idNow%2 == 1) etaNow *= coup2d;
1078     else              etaNow *= coup2u;   
1079     
1080     // Sum up contribution and return square of absolute value.
1081     eta += etaNow;
1082   }
1083   return (pow2(eta.real()) + pow2(eta.imag()));
1084
1085 }
1086
1087 //--------------------------------------------------------------------------
1088
1089 // Sum up quark, lepton, W+- and (for BSM) H+- loop contributions 
1090 // in Higgs -> gamma + gamma.
1091
1092 double ResonanceH::eta2gaga() {
1093
1094   // Initial values.
1095   complex eta = complex(0., 0.);
1096   int     idNow;
1097   double  ef, mLoop, epsilon, root, rootLog;
1098   complex phi, etaNow;
1099
1100   // Loop over s, c, b, t, mu, tau, W+-, H+- flavours.
1101   for (int idLoop = 0; idLoop < 8; ++idLoop) {
1102     if      (idLoop < 4) idNow = idLoop + 3;
1103     else if (idLoop < 6) idNow = 2 * idLoop + 5;
1104     else if (idLoop < 7) idNow = 24;
1105     else                 idNow = 37;
1106     if (idNow == 37 && higgsType == 0) continue;
1107  
1108     // Charge and loop integral parameter.
1109     ef      = (idNow < 20) ? couplingsPtr->ef(idNow) : 1.;
1110     mLoop   = (useRunLoopMass) ? particleDataPtr->mRun(idNow, mHat)
1111                                : particleDataPtr->m0(idNow);
1112     epsilon = pow2(2. * mLoop / mHat);
1113
1114     // Value of loop integral.
1115     if (epsilon <= 1.) {
1116       root    = sqrt(1. - epsilon);
1117       rootLog = (epsilon < 1e-4) ? log(4. / epsilon - 2.)
1118                 : log( (1. + root) / (1. - root) );
1119       phi     = complex( -0.25 * (pow2(rootLog) - pow2(M_PI)), 
1120                 0.5 * M_PI * rootLog );
1121     } 
1122     else phi  = complex( pow2( asin(1. / sqrt(epsilon)) ), 0.);
1123
1124     // Expressions for quarks and leptons that depend on Higgs type.
1125     if (idNow < 17) { 
1126       if (higgsType < 3) etaNow = -0.5 * epsilon 
1127         * (complex(1., 0.) + (1. - epsilon) * phi);
1128       else etaNow = -0.5 * epsilon * phi;    
1129       if (idNow < 7 && idNow%2 == 1) etaNow *= 3. * pow2(ef) * coup2d;
1130       else if (idNow < 7 )           etaNow *= 3. * pow2(ef) * coup2u;
1131       else                           etaNow *=      pow2(ef) * coup2l;
1132     } 
1133
1134     // Expression for W+-.
1135     else if (idNow == 24) etaNow = (complex(0.5 + 0.75 * epsilon, 0.)
1136       + 0.75 * epsilon * (2. - epsilon) * phi) * coup2W;  
1137  
1138     // Expression for H+-.
1139    else etaNow = (complex(epsilon, 0.) - epsilon * epsilon * phi)
1140      * pow2(mW / mHchg) * coup2Hchg;      
1141     
1142     // Sum up contribution and return square of absolute value.
1143     eta       += etaNow;
1144   }
1145   return (pow2(eta.real()) + pow2(eta.imag()));
1146
1147 }
1148
1149 //--------------------------------------------------------------------------
1150
1151 // Sum up quark, lepton, W+- and (for BSM) H+- loop contributions 
1152 // in Higgs -> gamma + Z0.
1153
1154 double ResonanceH::eta2gaZ() {
1155
1156   // Initial values.
1157   complex eta = complex(0., 0.);
1158   int     idNow;
1159   double  ef, vf, mLoop, epsilon, epsPrime, root, rootLog, asinEps;
1160   complex phi, psi, phiPrime, psiPrime, fXY, f1, etaNow;
1161
1162   // Loop over s, c, b, t, mu , tau, W+-, H+- flavours.
1163   for (int idLoop = 0; idLoop < 7; ++idLoop) {
1164     if      (idLoop < 4) idNow = idLoop + 3;
1165     else if (idLoop < 6) idNow = 2 * idLoop + 5;
1166     else if (idLoop < 7) idNow = 24;
1167     else                 idNow = 37;
1168
1169     // Electroweak charges and loop integral parameters.
1170     ef        = (idNow < 20) ? couplingsPtr->ef(idNow) : 1.;
1171     vf        = (idNow < 20) ? couplingsPtr->vf(idNow) : 0.;
1172     mLoop     = (useRunLoopMass) ? particleDataPtr->mRun(idNow, mHat)
1173                                  : particleDataPtr->m0(idNow);
1174     epsilon   = pow2(2. * mLoop / mHat);
1175     epsPrime  = pow2(2. * mLoop / mZ);
1176
1177     // Value of loop integral for epsilon = 4 m^2 / sHat.
1178     if (epsilon <= 1.) {
1179       root    = sqrt(1. - epsilon);
1180       rootLog = (epsilon < 1e-4) ? log(4. / epsilon - 2.)
1181                 : log( (1. + root) / (1. - root) );
1182       phi     = complex( -0.25 * (pow2(rootLog) - pow2(M_PI)), 
1183                 0.5 * M_PI * rootLog );
1184       psi     = 0.5 * root * complex( rootLog, -M_PI); 
1185     } else {
1186       asinEps = asin(1. / sqrt(epsilon));
1187       phi     = complex( pow2(asinEps), 0.);
1188       psi     = complex( sqrt(epsilon - 1.) * asinEps, 0.);
1189     }
1190
1191     // Value of loop integral for epsilonPrime = 4 m^2 / m_Z^2.
1192     if (epsPrime <= 1.) {
1193       root     = sqrt(1. - epsPrime);
1194       rootLog  = (epsPrime < 1e-4) ? log(4. / epsPrime - 2.)
1195                  : log( (1. + root) / (1. - root) );
1196       phiPrime = complex( -0.25 * (pow2(rootLog) - pow2(M_PI)), 
1197                           0.5 * M_PI * rootLog );
1198       psiPrime = 0.5 * root * complex( rootLog, -M_PI); 
1199     } else {
1200       asinEps  = asin(1. / sqrt(epsPrime));
1201       phiPrime = complex( pow2(asinEps), 0.);
1202       psiPrime = complex( sqrt(epsPrime - 1.) * asinEps, 0.);
1203     }
1204
1205     // Combine the two loop integrals.
1206     fXY = (epsilon * epsPrime / (8. * pow2(epsilon - epsPrime)))
1207       * ( complex(epsilon - epsPrime, 0) 
1208       + epsilon * epsPrime * (phi - phiPrime)
1209       + 2. * epsilon * (psi - psiPrime) );
1210     f1 = - (epsilon * epsPrime / (2. * (epsilon - epsPrime)))
1211       * (phi - phiPrime);    
1212
1213     // Expressions for quarks and leptons that depend on Higgs type.
1214     if (idNow < 17) { 
1215       etaNow = (higgsType < 3) ? -fXY + 0.25 * f1 : 0.25 * f1;
1216       if (idNow < 7 && idNow%2 == 1) etaNow *= 3. * ef * vf * coup2d;
1217       else if (idNow < 7)         etaNow *= 3. * ef * vf * coup2u;
1218       else                     etaNow *=      ef * vf * coup2l;
1219
1220     // Expression for W+-.
1221     } else if (idNow == 24) {
1222       double coef1  = 3. - sin2tW / cos2tW;
1223       double coefXY = (1. + 2. / epsilon) * sin2tW / cos2tW 
1224         - (5. + 2. / epsilon);
1225       etaNow = -cos2tW * (coef1 * f1 + coefXY * fXY) * coup2W; 
1226
1227     // Expression for H+-.
1228     } else etaNow = (1. - 2. * sin2tW) * fXY * pow2(mW / mHchg) 
1229       * coup2Hchg;     
1230     
1231     // Sum up contribution and return square of absolute value.
1232     eta += etaNow;
1233   }
1234   return ( (pow2(eta.real()) + pow2(eta.imag())) / (sin2tW * cos2tW) );
1235
1236 }
1237  
1238 //==========================================================================
1239
1240 // The ResonanceHchg class.
1241 // Derived class for H+- properties.
1242
1243 //--------------------------------------------------------------------------
1244
1245 // Initialize constants.
1246
1247 void ResonanceHchg::initConstants() {
1248
1249   // Locally stored properties and couplings.
1250   useCubicWidth = settingsPtr->flag("Higgs:cubicWidth");
1251   thetaWRat     = 1. / (8. * couplingsPtr->sin2thetaW());
1252   mW            = particleDataPtr->m0(24);
1253   tanBeta       = settingsPtr->parm("HiggsHchg:tanBeta");
1254   tan2Beta      = tanBeta * tanBeta;
1255   coup2H1W      = settingsPtr->parm("HiggsHchg:coup2H1W");
1256
1257 }
1258
1259 //--------------------------------------------------------------------------
1260
1261 // Calculate various common prefactors for the current mass.
1262
1263 void ResonanceHchg::calcPreFac(bool) {
1264
1265   // Common coupling factors.
1266   alpEM         = couplingsPtr->alphaEM(mHat * mHat);
1267   alpS          = couplingsPtr->alphaS(mHat * mHat);
1268   colQ          = 3. * (1. + alpS / M_PI);
1269   preFac        = alpEM * thetaWRat * pow3(mHat) / pow2(mW); 
1270
1271 }
1272
1273 //--------------------------------------------------------------------------
1274
1275 // Calculate width for currently considered channel.
1276
1277 void ResonanceHchg::calcWidth(bool) {
1278
1279   // Check that above threshold.
1280   if (ps == 0.) return;
1281
1282   // H+- decay to fermions involves running masses.
1283   if (id1Abs < 17 && (id1Abs < 7 || id1Abs > 10)) {
1284     double mRun1   = particleDataPtr->mRun(id1Abs, mHat);
1285     double mRun2   = particleDataPtr->mRun(id2Abs, mHat);
1286     double mrRunDn = pow2(mRun1 / mHat);
1287     double mrRunUp = pow2(mRun2 / mHat);
1288     if (id1Abs%2 == 0) swap( mrRunDn, mrRunUp);
1289
1290     // Width to fermions: couplings, kinematics, colour factor. 
1291     widNow = preFac * max( 0., (mrRunDn * tan2Beta + mrRunUp / tan2Beta) 
1292            * (1. - mrRunDn - mrRunUp) - 4. *mrRunDn * mrRunUp ) * ps;
1293     if (id1Abs < 7) widNow *= colQ;
1294   }
1295
1296   // H+- decay to h0 + W+-.
1297   else if (id1Abs == 25 && id2Abs == 24) 
1298     widNow    = 0.5 * preFac * pow3(ps) * pow2(coup2H1W);
1299
1300 }
1301  
1302 //==========================================================================
1303
1304 // The ResonanceZprime class.
1305 // Derived class for gamma*/Z0/Z'^0 properties.
1306
1307 //--------------------------------------------------------------------------
1308
1309 // Initialize constants.
1310
1311 void ResonanceZprime::initConstants() {
1312
1313   // Locally stored properties and couplings.
1314   gmZmode     = settingsPtr->mode("Zprime:gmZmode");
1315   sin2tW      = couplingsPtr->sin2thetaW();
1316   cos2tW      = 1. - sin2tW;
1317   thetaWRat   = 1. / (16. * sin2tW * cos2tW);
1318
1319   // Properties of Z resonance.
1320   mZ          = particleDataPtr->m0(23);
1321   GammaZ      = particleDataPtr->mWidth(23);
1322   m2Z         = mZ*mZ;
1323   GamMRatZ    = GammaZ / mZ;
1324
1325   // Ensure that arrays initially empty.
1326   for (int i = 0; i < 20; ++i) afZp[i] = 0.;
1327   for (int i = 0; i < 20; ++i) vfZp[i] = 0.;
1328   
1329   // Store first-generation axial and vector couplings.
1330   afZp[1]     = settingsPtr->parm("Zprime:ad");
1331   afZp[2]     = settingsPtr->parm("Zprime:au");
1332   afZp[11]    = settingsPtr->parm("Zprime:ae");
1333   afZp[12]    = settingsPtr->parm("Zprime:anue");
1334   vfZp[1]     = settingsPtr->parm("Zprime:vd");
1335   vfZp[2]     = settingsPtr->parm("Zprime:vu");
1336   vfZp[11]    = settingsPtr->parm("Zprime:ve");
1337   vfZp[12]    = settingsPtr->parm("Zprime:vnue");
1338
1339   // Second and third generation could be carbon copy of this...
1340   if (settingsPtr->flag("Zprime:universality")) {
1341     for (int i = 3; i <= 6; ++i) {
1342       afZp[i]    = afZp[i-2];
1343       vfZp[i]    = vfZp[i-2];
1344       afZp[i+10] = afZp[i+8];
1345       vfZp[i+10] = vfZp[i+8];
1346     }
1347
1348   // ... or could have different couplings.   
1349   } else {
1350     afZp[3]   = settingsPtr->parm("Zprime:as");
1351     afZp[4]   = settingsPtr->parm("Zprime:ac");
1352     afZp[5]   = settingsPtr->parm("Zprime:ab");
1353     afZp[6]   = settingsPtr->parm("Zprime:at");
1354     afZp[13]  = settingsPtr->parm("Zprime:amu");
1355     afZp[14]  = settingsPtr->parm("Zprime:anumu");
1356     afZp[15]  = settingsPtr->parm("Zprime:atau");
1357     afZp[16]  = settingsPtr->parm("Zprime:anutau");
1358     vfZp[3]   = settingsPtr->parm("Zprime:vs");
1359     vfZp[4]   = settingsPtr->parm("Zprime:vc");
1360     vfZp[5]   = settingsPtr->parm("Zprime:vb");
1361     vfZp[6]   = settingsPtr->parm("Zprime:vt");
1362     vfZp[13]  = settingsPtr->parm("Zprime:vmu");
1363     vfZp[14]  = settingsPtr->parm("Zprime:vnumu");
1364     vfZp[15]  = settingsPtr->parm("Zprime:vtau");
1365     vfZp[16]  = settingsPtr->parm("Zprime:vnutau");
1366   }
1367
1368   // Coupling for Z' -> W+ W-.
1369   coupZpWW    = settingsPtr->parm("Zprime:coup2WW");
1370
1371 }
1372
1373 //--------------------------------------------------------------------------
1374
1375 // Calculate various common prefactors for the current mass.
1376
1377 void ResonanceZprime::calcPreFac(bool calledFromInit) {
1378
1379   // Common coupling factors.
1380   alpEM       = couplingsPtr->alphaEM(mHat * mHat);
1381   alpS        = couplingsPtr->alphaS(mHat * mHat);
1382   colQ        = 3. * (1. + alpS / M_PI);
1383   preFac      = alpEM * thetaWRat * mHat / 3.;
1384
1385   // When call for incoming flavour need to consider gamma*/Z0 mix.
1386   if (!calledFromInit) {
1387
1388     // Couplings when an incoming fermion is specified; elso only pure Z'0.
1389     ei2       = 0.;
1390     eivi      = 0.;  
1391     vai2      = 0.;
1392     eivpi     = 0.;
1393     vaivapi   = 0.,
1394     vapi2     = 1.;
1395     int idInFlavAbs = abs(idInFlav);
1396     if (idInFlavAbs > 0 && idInFlavAbs < 19) {
1397       double ei  = couplingsPtr->ef(idInFlavAbs);
1398       double ai  = couplingsPtr->af(idInFlavAbs);
1399       double vi  = couplingsPtr->vf(idInFlavAbs);
1400       double api = afZp[idInFlavAbs];
1401       double vpi = vfZp[idInFlavAbs];
1402       ei2     = ei * ei;
1403       eivi    = ei * vi;
1404       vai2    = vi * vi + ai * ai; 
1405       eivpi   = ei * vpi;
1406       vaivapi = vi * vpi + ai * api;; 
1407       vapi2   = vpi * vpi + api * api;
1408     }
1409
1410     // Calculate prefactors for gamma/interference/Z0 terms.
1411     double sH     = mHat * mHat;
1412     double propZ  = sH / ( pow2(sH - m2Z) + pow2(sH * GamMRatZ) );
1413     double propZp = sH / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1414     gamNorm   = ei2;
1415     gamZNorm  = 2. * eivi * thetaWRat * (sH - m2Z) * propZ;
1416     ZNorm     = vai2 * pow2(thetaWRat) * sH * propZ;
1417     gamZpNorm = 2. * eivpi * thetaWRat * (sH - m2Res) * propZp;
1418     ZZpNorm   = 2. * vaivapi * pow2(thetaWRat) * ((sH - m2Res) * (sH - m2Z) 
1419               + sH * GamMRat * sH * GamMRatZ) * propZ * propZp;
1420     ZpNorm    = vapi2 * pow2(thetaWRat) * sH * propZp;
1421
1422     // Rescale Z0 height normalization to compensate for a width one??
1423     //if (doForceWidth) {
1424     //  intNorm *= forceFactor;
1425     //  resNorm *= forceFactor;
1426     //}
1427
1428     // Optionally only keep some of gamma*, Z0 and Z' terms.
1429     if (gmZmode == 1) {gamZNorm = 0; ZNorm = 0.; gamZpNorm = 0.; 
1430       ZZpNorm = 0.; ZpNorm = 0.;}
1431     if (gmZmode == 2) {gamNorm = 0.; gamZNorm = 0.; gamZpNorm = 0.; 
1432       ZZpNorm = 0.; ZpNorm = 0.;}
1433     if (gmZmode == 3) {gamNorm = 0.; gamZNorm = 0.; ZNorm = 0.;
1434       gamZpNorm = 0.; ZZpNorm = 0.;}
1435     if (gmZmode == 4) {gamZpNorm = 0.; ZZpNorm = 0.; ZpNorm = 0.;}
1436     if (gmZmode == 5) {gamZNorm = 0.; ZNorm = 0.; ZZpNorm = 0.;}
1437     if (gmZmode == 6) {gamNorm = 0.; gamZNorm = 0.; gamZpNorm = 0.;}
1438   }
1439
1440 }
1441
1442 //--------------------------------------------------------------------------
1443
1444 // Calculate width for currently considered channel.
1445
1446 void ResonanceZprime::calcWidth(bool calledFromInit) {
1447
1448   // Check that above threshold. 
1449   if (ps == 0.) return;
1450   
1451   // At initialization only the pure Z'0 should be considered.
1452   if (calledFromInit) {
1453
1454     // Contributions from three fermion generations.
1455     if ( id1Abs < 7 || (id1Abs > 10 && id1Abs < 17) ) {
1456       double apf = afZp[id1Abs];
1457       double vpf = vfZp[id1Abs];
1458       widNow = preFac * ps * (vpf*vpf * (1. + 2. * mr1) 
1459              + apf*apf * ps*ps); 
1460       if (id1Abs < 7) widNow *= colQ;
1461
1462     // Contribution from Z'0 -> W^+ W^-.
1463     } else if (id1Abs == 24) {
1464       widNow = preFac * pow2(coupZpWW * cos2tW) * pow3(ps)
1465         * (1. + mr1*mr1 + mr2*mr2 + 10. * (mr1 + mr2 + mr1 * mr2));
1466     }
1467   }   
1468
1469   // When call for incoming flavour need to consider full mix.
1470   else {
1471
1472     // Contributions from three fermion generations.
1473     if ( id1Abs < 7 || (id1Abs > 10 && id1Abs < 17) ) {
1474
1475       // Couplings of gamma^*/Z^0/Z'^0  to final flavour 
1476       double ef  = couplingsPtr->ef(id1Abs);
1477       double af  = couplingsPtr->af(id1Abs);
1478       double vf  = couplingsPtr->vf(id1Abs);
1479       double apf = afZp[id1Abs];
1480       double vpf = vfZp[id1Abs];
1481
1482       // Combine couplings with kinematical factors.
1483       double kinFacA  = pow3(ps);
1484       double kinFacV  = ps * (1. + 2. * mr1);
1485       double ef2      = ef * ef * kinFacV;
1486       double efvf     = ef * vf * kinFacV;
1487       double vaf2     = vf * vf * kinFacV + af * af * kinFacA; 
1488       double efvpf    = ef * vpf * kinFacV;
1489       double vafvapf  = vf * vpf * kinFacV + af * apf * kinFacA;
1490       double vapf2    = vpf * vpf * kinFacV + apf * apf * kinFacA;
1491     
1492       // Relative outwidths: combine instate, propagator and outstate.
1493       widNow = gamNorm * ef2 + gamZNorm * efvf + ZNorm * vaf2
1494              + gamZpNorm * efvpf + ZZpNorm * vafvapf + ZpNorm * vapf2;
1495       if (id1Abs < 7) widNow *= colQ;
1496
1497     // Contribution from Z'0 -> W^+ W^-.
1498     } else if (id1Abs == 24) {
1499       widNow = ZpNorm * pow2(coupZpWW * cos2tW) * pow3(ps)
1500         * (1. + mr1*mr1 + mr2*mr2 + 10. * (mr1 + mr2 + mr1 * mr2));
1501     }
1502   }
1503
1504 }
1505  
1506 //==========================================================================
1507
1508 // The ResonanceWprime class.
1509 // Derived class for W'+- properties.
1510
1511 //--------------------------------------------------------------------------
1512
1513 // Initialize constants.
1514
1515 void ResonanceWprime::initConstants() {
1516
1517   // Locally stored properties and couplings.
1518   thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());
1519   cos2tW    = couplingsPtr->cos2thetaW();
1520
1521   // Axial and vector couplings of fermions.
1522   aqWp      = settingsPtr->parm("Wprime:aq");
1523   vqWp      = settingsPtr->parm("Wprime:vq");
1524   alWp      = settingsPtr->parm("Wprime:al");
1525   vlWp      = settingsPtr->parm("Wprime:vl");
1526
1527   // Coupling for W' -> W Z.
1528   coupWpWZ    = settingsPtr->parm("Wprime:coup2WZ");
1529
1530 }
1531
1532 //--------------------------------------------------------------------------
1533
1534 // Calculate various common prefactors for the current mass.
1535
1536 void ResonanceWprime::calcPreFac(bool) {
1537
1538   // Common coupling factors.
1539   alpEM     = couplingsPtr->alphaEM(mHat * mHat);
1540   alpS      = couplingsPtr->alphaS(mHat * mHat);
1541   colQ      = 3. * (1. + alpS / M_PI);
1542   preFac    = alpEM * thetaWRat * mHat;
1543
1544 }
1545
1546 //--------------------------------------------------------------------------
1547
1548 // Calculate width for currently considered channel.
1549
1550 void ResonanceWprime::calcWidth(bool) {
1551
1552   // Check that above threshold.
1553   if (ps == 0.) return;
1554
1555   // Decay to quarks involves colour factor and CKM matrix.
1556   if (id1Abs > 0 && id1Abs < 7) widNow    
1557     = preFac * ps * 0.5 * ((vqWp*vqWp + aqWp * aqWp) 
1558     + 6. * (vqWp*vqWp - aqWp * aqWp) * sqrt(mr1 *mr2)) 
1559     * (1. - 0.5 * (mr1 + mr2) - 0.5 * pow2(mr1 - mr2))
1560     * colQ * couplingsPtr->V2CKMid(id1Abs, id2Abs);
1561
1562   // Decay to leptons simpler. 
1563   else if (id1Abs > 10 && id1Abs < 17) widNow    
1564     = preFac * ps * 0.5 * ((vlWp*vqWp + alWp * aqWp) 
1565     + 6. * (vlWp*vqWp - alWp * aqWp) * sqrt(mr1 *mr2))  
1566     * (1. - 0.5 * (mr1 + mr2) - 0.5 * pow2(mr1 - mr2));
1567
1568   // Decay to W^+- Z^0.
1569   else if (id1Abs == 24 && id2Abs == 23) widNow
1570     = preFac * 0.25 * pow2(coupWpWZ) * cos2tW * (mr1 / mr2) * pow3(ps)
1571     * (1. + mr1*mr1 + mr2*mr2 + 10. * (mr1 + mr2 + mr1 * mr2));
1572
1573 }
1574  
1575 //==========================================================================
1576
1577 // The ResonanceRhorizontal class.
1578 // Derived class for R^0 (horizontal gauge boson) properties.
1579
1580 //--------------------------------------------------------------------------
1581
1582 // Initialize constants.
1583
1584 void ResonanceRhorizontal::initConstants() {
1585
1586   // Locally stored properties and couplings.
1587   thetaWRat = 1. / (12. * couplingsPtr->sin2thetaW());
1588
1589 }
1590
1591 //--------------------------------------------------------------------------
1592
1593 // Calculate various common prefactors for the current mass.
1594
1595 void ResonanceRhorizontal::calcPreFac(bool) {
1596
1597   // Common coupling factors.
1598   alpEM     = couplingsPtr->alphaEM(mHat * mHat);
1599   alpS      = couplingsPtr->alphaS(mHat * mHat);
1600   colQ      = 3. * (1. + alpS / M_PI);
1601   preFac    = alpEM * thetaWRat * mHat;
1602
1603 }
1604
1605 //--------------------------------------------------------------------------
1606
1607 // Calculate width for currently considered channel.
1608
1609 void ResonanceRhorizontal::calcWidth(bool) {
1610
1611   // Check that above threshold.
1612   if (ps == 0.) return;
1613
1614   // R -> f fbar. Colour factor for quarks.
1615   widNow    = preFac * ps * (2. - mr1 - mr2 - pow2(mr1 - mr2));
1616   if (id1Abs < 9) widNow *= colQ;
1617
1618 }
1619  
1620 //==========================================================================
1621
1622 // The ResonanceExcited class.
1623 // Derived class for excited-fermion properties.
1624
1625 //--------------------------------------------------------------------------
1626
1627 // Initialize constants.
1628
1629 void ResonanceExcited::initConstants() {
1630
1631   // Locally stored properties and couplings.
1632   Lambda        = settingsPtr->parm("ExcitedFermion:Lambda");
1633   coupF         = settingsPtr->parm("ExcitedFermion:coupF");
1634   coupFprime    = settingsPtr->parm("ExcitedFermion:coupFprime");
1635   coupFcol      = settingsPtr->parm("ExcitedFermion:coupFcol");
1636   sin2tW        = couplingsPtr->sin2thetaW();
1637   cos2tW        = 1. - sin2tW;
1638
1639 }
1640
1641 //--------------------------------------------------------------------------
1642
1643 // Calculate various common prefactors for the current mass.
1644
1645 void ResonanceExcited::calcPreFac(bool) {
1646
1647   // Common coupling factors.
1648   alpEM         = couplingsPtr->alphaEM(mHat * mHat);
1649   alpS          = couplingsPtr->alphaS(mHat * mHat);
1650   preFac        = pow3(mHat) / pow2(Lambda);
1651
1652 }
1653
1654 //--------------------------------------------------------------------------
1655
1656 // Calculate width for currently considered channel.
1657
1658 void ResonanceExcited::calcWidth(bool) {
1659
1660   // Check that above threshold.
1661   if (ps == 0.) return;
1662
1663   // f^* -> f g.
1664   if (id1Abs == 21) widNow = preFac * alpS * pow2(coupFcol) / 3.;
1665
1666   // f^* -> f gamma.
1667   else if (id1Abs == 22) {
1668     double chgI3 = (id2Abs%2 == 0) ? 0.5 : -0.5;
1669     double chgY  = (id2Abs < 9) ? 1. / 6. : -0.5; 
1670     double chg   = chgI3 * coupF + chgY * coupFprime;
1671     widNow       = preFac * alpEM * pow2(chg) / 4.;      
1672   }
1673
1674   // f^* -> f Z^0.
1675   else if (id1Abs == 23) {
1676     double chgI3 = (id2Abs%2 == 0) ? 0.5 : -0.5;
1677     double chgY  = (id2Abs < 9) ? 1. / 6. : -0.5; 
1678     double chg   = chgI3 * cos2tW * coupF - chgY * sin2tW * coupFprime;
1679     widNow       = preFac * (alpEM * pow2(chg) / (8. * sin2tW * cos2tW))
1680                  * ps*ps * (2. + mr1);      
1681   }
1682
1683   // f^* -> f' W^+-.
1684   else if (id1Abs == 24) widNow  = preFac * (alpEM * pow2(coupF) 
1685                  / (16. * sin2tW)) * ps*ps * (2. + mr1);      
1686
1687 }
1688  
1689 //==========================================================================
1690
1691 // The ResonanceGraviton class.
1692 // Derived class for excited Graviton properties.
1693
1694 //--------------------------------------------------------------------------
1695
1696 // Initialize constants.
1697
1698 void ResonanceGraviton::initConstants() {
1699
1700   // SMinBulk = off/on, use universal coupling (kappaMG) 
1701   // or individual (Gxx) between graviton and SM particles.
1702   eDsmbulk   = settingsPtr->flag("ExtraDimensionsG*:SMinBulk");
1703   eDvlvl = false;
1704   if (eDsmbulk) eDvlvl = settingsPtr->flag("ExtraDimensionsG*:VLVL");
1705   kappaMG    = settingsPtr->parm("ExtraDimensionsG*:kappaMG");
1706   for (int i = 0; i < 27; ++i) eDcoupling[i] = 0.;
1707   double tmp_coup = settingsPtr->parm("ExtraDimensionsG*:Gqq");
1708   for (int i = 1; i <= 4; ++i)  eDcoupling[i] = tmp_coup;
1709   eDcoupling[5] = settingsPtr->parm("ExtraDimensionsG*:Gbb"); 
1710   eDcoupling[6] = settingsPtr->parm("ExtraDimensionsG*:Gtt");
1711   tmp_coup = settingsPtr->parm("ExtraDimensionsG*:Gll");
1712   for (int i = 11; i <= 16; ++i) eDcoupling[i] = tmp_coup;
1713   eDcoupling[21] = settingsPtr->parm("ExtraDimensionsG*:Ggg");
1714   eDcoupling[22] = settingsPtr->parm("ExtraDimensionsG*:Ggmgm");
1715   eDcoupling[23] = settingsPtr->parm("ExtraDimensionsG*:GZZ");
1716   eDcoupling[24] = settingsPtr->parm("ExtraDimensionsG*:GWW");
1717   eDcoupling[25] = settingsPtr->parm("ExtraDimensionsG*:Ghh");
1718
1719 }
1720
1721 //--------------------------------------------------------------------------
1722
1723 // Calculate various common prefactors for the current mass.
1724
1725 void ResonanceGraviton::calcPreFac(bool) {
1726
1727   // Common coupling factors.
1728   alpS          = couplingsPtr->alphaS(mHat * mHat);
1729   colQ          = 3. * (1. + alpS / M_PI);
1730   preFac        = mHat / M_PI;
1731
1732 }
1733
1734 //--------------------------------------------------------------------------
1735
1736 // Calculate width for currently considered channel.
1737
1738 void ResonanceGraviton::calcWidth(bool) {
1739
1740   // Check that above threshold.
1741   if (ps == 0.) return;
1742
1743   // Widths to fermion pairs.
1744   if (id1Abs < 19) {
1745      widNow  = preFac * pow3(ps) * (1. + 8. * mr1 / 3.) / 320.;        
1746      if (id1Abs < 9) widNow *= colQ;
1747       
1748   // Widths to gluon and photon pair.
1749   } else if (id1Abs == 21) {
1750     widNow = preFac / 20.;
1751   } else if (id1Abs == 22) {
1752     widNow = preFac / 160.;
1753      
1754   // Widths to Z0 Z0 and W+ W- pair.
1755   } else if (id1Abs == 23 || id1Abs == 24) {
1756     // Longitudinal W/Z only.
1757     if (eDvlvl) {
1758       widNow = preFac * pow(ps,5) / 480.;
1759     // Transverse W/Z contributions as well.
1760     } else {
1761       widNow  = preFac * ps * (13. / 12. + 14. * mr1 / 3. + 4. * mr1 * mr1) 
1762               / 80.;
1763     }
1764     if (id1Abs == 23) widNow *= 0.5;
1765
1766   // Widths to h h pair.
1767   } else if (id1Abs == 25) {
1768     widNow = preFac * pow(ps,5) / 960.;
1769   }
1770
1771   // RS graviton coupling
1772   if (eDsmbulk) widNow *= 2. * pow2(eDcoupling[min( id1Abs, 26)] * mHat);  
1773   else          widNow *= pow2(kappaMG);
1774
1775 }
1776
1777 //==========================================================================
1778
1779 // The ResonanceKKgluon class.
1780 // Derived class for excited kk-gluon properties.
1781
1782 //--------------------------------------------------------------------------
1783
1784 // Initialize constants.
1785
1786 void ResonanceKKgluon::initConstants() {
1787
1788   // KK-gluon gv/ga couplings and interference.
1789   for (int i = 0; i < 10; ++i) { eDgv[i] = 0.; eDga[i] = 0.; }
1790   double tmp_gL = settingsPtr->parm("ExtraDimensionsG*:KKgqL");
1791   double tmp_gR = settingsPtr->parm("ExtraDimensionsG*:KKgqR");
1792   for (int i = 1; i <= 4; ++i) { 
1793     eDgv[i] = 0.5 * (tmp_gL + tmp_gR);
1794     eDga[i] = 0.5 * (tmp_gL - tmp_gR); 
1795   }
1796   tmp_gL = settingsPtr->parm("ExtraDimensionsG*:KKgbL"); 
1797   tmp_gR = settingsPtr->parm("ExtraDimensionsG*:KKgbR"); 
1798   eDgv[5] = 0.5 * (tmp_gL + tmp_gR); eDga[5] = 0.5 * (tmp_gL - tmp_gR); 
1799   tmp_gL = settingsPtr->parm("ExtraDimensionsG*:KKgtL"); 
1800   tmp_gR = settingsPtr->parm("ExtraDimensionsG*:KKgtR"); 
1801   eDgv[6] = 0.5 * (tmp_gL + tmp_gR); eDga[6] = 0.5 * (tmp_gL - tmp_gR); 
1802   interfMode    = settingsPtr->mode("ExtraDimensionsG*:KKintMode");
1803
1804 }
1805
1806 //--------------------------------------------------------------------------
1807
1808 // Calculate various common prefactors for the current mass.
1809
1810 void ResonanceKKgluon::calcPreFac(bool calledFromInit) {
1811
1812   // Common coupling factors.
1813   alpS          = couplingsPtr->alphaS(mHat * mHat);
1814   preFac        = alpS * mHat / 6;
1815
1816   // When call for incoming flavour need to consider g*/gKK mix.
1817   if (!calledFromInit) {
1818     // Calculate prefactors for g/interference/gKK terms.
1819     int idInFlavAbs = abs(idInFlav);
1820     double sH = mHat * mHat;
1821     normSM   = 1;
1822     normInt  = 2. * eDgv[min(idInFlavAbs, 9)] * sH * (sH - m2Res)
1823               / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1824     normKK   = ( pow2(eDgv[min(idInFlavAbs, 9)]) 
1825                + pow2(eDga[min(idInFlavAbs, 9)]) ) * sH * sH 
1826               / ( pow2(sH - m2Res) + pow2(sH * GamMRat) );
1827
1828     // Optionally only keep g* or gKK term.
1829     if (interfMode == 1) {normInt = 0.; normKK = 0.;}
1830     if (interfMode == 2) {normSM = 0.; normInt = 0.; normKK = 1.;}
1831   }
1832
1833 }
1834
1835 //--------------------------------------------------------------------------
1836
1837 // Calculate width for currently considered channel.
1838
1839 void ResonanceKKgluon::calcWidth(bool calledFromInit) {
1840
1841   // Check that above threshold.
1842   if (ps == 0.) return;
1843
1844   // Widths to quark pairs.
1845   if (id1Abs > 9) return;
1846
1847   if (calledFromInit) {
1848     widNow = preFac * ps * (pow2(eDgv[min(id1Abs, 9)]) * (1. + 2.*mr1) 
1849                          +  pow2(eDga[min(id1Abs, 9)]) * (1. - 4.*mr1) );
1850   } else {
1851     // Relative outwidths: combine instate, propagator and outstate.
1852     widNow = normSM  * ps * (1. + 2. * mr1) 
1853            + normInt * ps * eDgv[min(id1Abs, 9)] * (1. + 2. * mr1)
1854            + normKK  * ps * (pow2(eDgv[min(id1Abs, 9)]) * (1. + 2.*mr1) 
1855                           +  pow2(eDga[min(id1Abs, 9)]) * (1. - 4.*mr1) );
1856     widNow *= preFac;
1857   }
1858
1859
1860
1861 //==========================================================================
1862
1863 // The ResonanceLeptoquark class.
1864 // Derived class for leptoquark properties.
1865
1866 //--------------------------------------------------------------------------
1867
1868 // Initialize constants.
1869
1870 void ResonanceLeptoquark::initConstants() {
1871
1872   // Locally stored properties and couplings.
1873   kCoup      = settingsPtr->parm("LeptoQuark:kCoup");
1874
1875   // Check that flavour info in decay channel is correctly set.
1876   int id1Now = particlePtr->channel(0).product(0);
1877   int id2Now = particlePtr->channel(0).product(1);
1878   if (id1Now < 1 || id1Now > 5) {
1879     infoPtr->errorMsg("Error in ResonanceLeptoquark::init:"
1880       " unallowed input quark flavour reset to u"); 
1881     id1Now   = 2;
1882     particlePtr->channel(0).product(0, id1Now);
1883   }
1884   if (abs(id2Now) < 11 || abs(id2Now) > 16) {
1885     infoPtr->errorMsg("Error in ResonanceLeptoquark::init:"
1886       " unallowed input lepton flavour reset to e-"); 
1887     id2Now   = 11;
1888     particlePtr->channel(0).product(1, id2Now);
1889   }
1890
1891   // Set/overwrite charge and name of particle.
1892   bool changed  = particlePtr->hasChanged();
1893   int chargeLQ  = particleDataPtr->chargeType(id1Now) 
1894                 + particleDataPtr->chargeType(id2Now);
1895   particlePtr->setChargeType(chargeLQ); 
1896   string nameLQ = "LQ_" + particleDataPtr->name(id1Now) + ","
1897                 + particleDataPtr->name(id2Now);
1898   particlePtr->setNames(nameLQ, nameLQ + "bar"); 
1899   if (!changed) particlePtr->setHasChanged(false);
1900
1901 }
1902
1903 //--------------------------------------------------------------------------
1904
1905 // Calculate various common prefactors for the current mass.
1906
1907 void ResonanceLeptoquark::calcPreFac(bool) {
1908
1909   // Common coupling factors.
1910   alpEM         = couplingsPtr->alphaEM(mHat * mHat);
1911   preFac        = 0.25 * alpEM * kCoup * mHat; 
1912
1913 }
1914
1915 //--------------------------------------------------------------------------
1916
1917 // Calculate width for currently considered channel.
1918
1919 void ResonanceLeptoquark::calcWidth(bool) {
1920
1921   // Check that above threshold.
1922   if (ps == 0.) return;
1923   
1924   // Width into lepton plus quark.
1925   if (id1Abs > 10 && id1Abs < 17 && id2Abs < 7) widNow = preFac * pow3(ps);
1926
1927 }
1928  
1929 //==========================================================================
1930
1931 // The ResonanceNuRight class.
1932 // Derived class for righthanded Majorana neutrino properties.
1933
1934 //--------------------------------------------------------------------------
1935
1936 // Initialize constants.
1937
1938 void ResonanceNuRight::initConstants() {
1939
1940   // Locally stored properties and couplings: righthanded W mass.
1941   thetaWRat = 1. / (768. * M_PI * pow2(couplingsPtr->sin2thetaW()));
1942   mWR       = particleDataPtr->m0(9900024); 
1943
1944 }
1945
1946 //--------------------------------------------------------------------------
1947
1948 // Calculate various common prefactors for the current mass.
1949
1950 void ResonanceNuRight::calcPreFac(bool) {
1951
1952   // Common coupling factors.
1953   alpEM     = couplingsPtr->alphaEM(mHat * mHat);
1954   alpS      = couplingsPtr->alphaS(mHat * mHat);
1955   colQ      = 3. * (1. + alpS / M_PI);
1956   preFac    = pow2(alpEM) * thetaWRat * pow5(mHat) / pow4(max(mHat, mWR));
1957
1958 }
1959
1960 //--------------------------------------------------------------------------
1961
1962 // Calculate width for currently considered channel.
1963
1964 void ResonanceNuRight::calcWidth(bool) {
1965
1966   // Check that above threshold.
1967   if (mHat < mf1 + mf2 + mf3 + MASSMARGIN) return;
1968
1969   // Coupling part of widths to l- q qbar', l- l'+ nu_lR' and c.c.
1970   widNow    = (id2Abs < 9 && id3Abs < 9) 
1971             ? preFac * colQ * couplingsPtr->V2CKMid(id2, id3) : preFac;  
1972
1973   // Phase space corrections in decay. Must have y < 1.
1974   double x  = (mf1 + mf2 + mf3) / mHat; 
1975   double x2 = x * x;
1976   double fx = 1. - 8. * x2 + 8. * pow3(x2) - pow4(x2) 
1977             - 24. * pow2(x2) * log(x);
1978   double y  = min( 0.999, pow2(mHat / mWR) );
1979   double fy = ( 12. * (1. - y) * log(1. - y) + 12. * y - 6. * y*y 
1980             - 2.* pow3(y) ) / pow4(y);
1981   widNow   *= fx * fy;
1982
1983 }
1984  
1985 //==========================================================================
1986
1987 // The ResonanceZRight class.
1988 // Derived class for Z_R^0 properties.
1989
1990 //--------------------------------------------------------------------------
1991
1992 // Initialize constants.
1993
1994 void ResonanceZRight::initConstants() {
1995
1996   // Locally stored properties and couplings: righthanded W mass.
1997   sin2tW    = couplingsPtr->sin2thetaW();
1998   thetaWRat = 1. / (48. * sin2tW  * (1. - sin2tW) * (1. - 2. * sin2tW) );
1999
2000 }
2001
2002 //--------------------------------------------------------------------------
2003
2004 // Calculate various common prefactors for the current mass.
2005
2006 void ResonanceZRight::calcPreFac(bool) {
2007
2008   // Common coupling factors.
2009   alpEM     = couplingsPtr->alphaEM(mHat * mHat);
2010   alpS      = couplingsPtr->alphaS(mHat * mHat);
2011   colQ      = 3. * (1. + alpS / M_PI);
2012   preFac    = alpEM * thetaWRat * mHat;
2013
2014 }
2015
2016 //--------------------------------------------------------------------------
2017
2018 // Calculate width for currently considered channel.
2019
2020 void ResonanceZRight::calcWidth(bool) {
2021
2022   // Check that above threshold.
2023   if (ps == 0.) return;
2024
2025   // Couplings to q qbar and l+ l-.
2026   double vf = 0.;
2027   double af = 0.;
2028   double symMaj = 1.;
2029   if (id1Abs < 9 && id1Abs%2 == 1) {
2030     af      = -1. + 2. * sin2tW; 
2031     vf      = -1. + 4. * sin2tW / 3.;
2032   } else if (id1Abs < 9) {   
2033     af      = 1. - 2. * sin2tW; 
2034     vf      = 1. - 8. * sin2tW / 3.;
2035   } else if (id1Abs < 19 && id1Abs%2 == 1) {  
2036     af      = -1. + 2. * sin2tW; 
2037     vf      = -1. + 4. * sin2tW;
2038        
2039   // Couplings to nu_L nu_Lbar and nu_R nu_Rbar, both assumed Majoranas.
2040   } else if (id1Abs < 19) {
2041     af      = -2. * sin2tW; 
2042     vf      = 0.;
2043     symMaj  = 0.5; 
2044   } else {
2045     af      = 2. * (1. - sin2tW); 
2046     vf      = 0.;
2047     symMaj  = 0.5; 
2048   }
2049
2050   // Width expression, including phase space and colour factor.
2051   widNow    = preFac * (vf*vf * (1. + 2. * mr1) + af*af * ps*ps) * ps 
2052             * symMaj;
2053   if (id1Abs < 9) widNow *= colQ;
2054
2055 }
2056  
2057 //==========================================================================
2058
2059 // The ResonanceWRight class.
2060 // Derived class for W_R+- properties.
2061
2062 //--------------------------------------------------------------------------
2063
2064 // Initialize constants.
2065
2066 void ResonanceWRight::initConstants() {
2067
2068   // Locally stored properties and couplings.
2069   thetaWRat     = 1. / (12. * couplingsPtr->sin2thetaW());
2070
2071 }
2072
2073 //--------------------------------------------------------------------------
2074
2075 // Calculate various common prefactors for the current mass.
2076
2077 void ResonanceWRight::calcPreFac(bool) {
2078
2079   // Common coupling factors.
2080   alpEM     = couplingsPtr->alphaEM(mHat * mHat);
2081   alpS      = couplingsPtr->alphaS(mHat * mHat);
2082   colQ      = 3. * (1. + alpS / M_PI);
2083   preFac    = alpEM * thetaWRat * mHat;
2084
2085 }
2086
2087 //--------------------------------------------------------------------------
2088
2089 // Calculate width for currently considered channel.
2090
2091 void ResonanceWRight::calcWidth(bool) {
2092
2093   // Check that above threshold.
2094   if (ps == 0.) return;
2095
2096   // Combine kinematics with colour factor and CKM couplings.
2097   widNow    = preFac * (1. - 0.5 * (mr1 + mr2) - 0.5 * pow2(mr1 - mr2))
2098             * ps;
2099   if (id1Abs < 9) widNow *= colQ * couplingsPtr->V2CKMid(id1Abs, id2Abs);
2100
2101 }
2102  
2103 //==========================================================================
2104
2105 // The ResonanceHchgchgLeft class.
2106 // Derived class for H++/H-- (left) properties.
2107
2108 //--------------------------------------------------------------------------
2109
2110 // Initialize constants.
2111
2112 void ResonanceHchgchgLeft::initConstants() {
2113
2114   // Read in Yukawa matrix for couplings to a lepton pair.
2115   yukawa[1][1]  = settingsPtr->parm("LeftRightSymmmetry:coupHee");
2116   yukawa[2][1]  = settingsPtr->parm("LeftRightSymmmetry:coupHmue");
2117   yukawa[2][2]  = settingsPtr->parm("LeftRightSymmmetry:coupHmumu");
2118   yukawa[3][1]  = settingsPtr->parm("LeftRightSymmmetry:coupHtaue");
2119   yukawa[3][2]  = settingsPtr->parm("LeftRightSymmmetry:coupHtaumu");
2120   yukawa[3][3]  = settingsPtr->parm("LeftRightSymmmetry:coupHtautau");
2121   
2122   // Locally stored properties and couplings.
2123   gL            = settingsPtr->parm("LeftRightSymmmetry:gL");
2124   vL            = settingsPtr->parm("LeftRightSymmmetry:vL");
2125   mW            = particleDataPtr->m0(24);
2126
2127 }
2128
2129 //--------------------------------------------------------------------------
2130
2131 // Calculate various common prefactors for the current mass.
2132
2133 void ResonanceHchgchgLeft::calcPreFac(bool) {
2134
2135   // Common coupling factors.
2136   preFac        = mHat / (8. * M_PI);
2137
2138 }
2139
2140 //--------------------------------------------------------------------------
2141
2142 // Calculate width for currently considered channel.
2143
2144 void ResonanceHchgchgLeft::calcWidth(bool) {
2145
2146   // Check that above threshold.
2147   if (ps == 0.) return;
2148
2149   // H++-- width to a pair of leptons. Combinatorial factor of 2.
2150   if (id1Abs < 17 && id2Abs < 17) {
2151     widNow    = preFac * pow2(yukawa[(id1Abs-9)/2][(id2Abs-9)/2]) * ps;
2152     if (id2Abs != id1Abs) widNow *= 2.;
2153   }
2154       
2155   // H++-- width to a pair of lefthanded W's.
2156   else if (id1Abs == 24 && id2Abs == 24)  
2157     widNow    = preFac * 0.5 * pow2(gL*gL * vL / mW) 
2158               * (3. * mr1 + 0.25 / mr1 - 1.) * ps;
2159
2160 }
2161  
2162 //==========================================================================
2163
2164 // The ResonanceHchgchgRight class.
2165 // Derived class for H++/H-- (right) properties.
2166
2167 //--------------------------------------------------------------------------
2168
2169 // Initialize constants.
2170
2171 void ResonanceHchgchgRight::initConstants() {
2172
2173   // Read in Yukawa matrix for couplings to a lepton pair.
2174   yukawa[1][1]  = settingsPtr->parm("LeftRightSymmmetry:coupHee");
2175   yukawa[2][1]  = settingsPtr->parm("LeftRightSymmmetry:coupHmue");
2176   yukawa[2][2]  = settingsPtr->parm("LeftRightSymmmetry:coupHmumu");
2177   yukawa[3][1]  = settingsPtr->parm("LeftRightSymmmetry:coupHtaue");
2178   yukawa[3][2]  = settingsPtr->parm("LeftRightSymmmetry:coupHtaumu");
2179   yukawa[3][3]  = settingsPtr->parm("LeftRightSymmmetry:coupHtautau");
2180   
2181   // Locally stored properties and couplings.
2182   idWR          = 9000024;
2183   gR            = settingsPtr->parm("LeftRightSymmmetry:gR");
2184
2185 }
2186
2187 //--------------------------------------------------------------------------
2188
2189 // Calculate various common prefactors for the current mass.
2190
2191 void ResonanceHchgchgRight::calcPreFac(bool) {
2192
2193   // Common coupling factors.
2194   preFac        = mHat / (8. * M_PI);
2195
2196 }
2197
2198 //--------------------------------------------------------------------------
2199
2200 // Calculate width for currently considered channel.
2201
2202 void ResonanceHchgchgRight::calcWidth(bool) {
2203
2204   // Check that above threshold.
2205   if (ps == 0.) return;
2206
2207   // H++-- width to a pair of leptons. Combinatorial factor of 2.
2208   if (id1Abs < 17 && id2Abs < 17) {
2209     widNow    = preFac * pow2(yukawa[(id1Abs-9)/2][(id2Abs-9)/2]) * ps;
2210     if (id2Abs != id1Abs) widNow *= 2.;
2211   }
2212       
2213   // H++-- width to a pair of lefthanded W's.
2214   else if (id1Abs == idWR && id2Abs == idWR)  
2215     widNow    = preFac * pow2(yukawa[(id1Abs-9)/2][(id2Abs-9)/2]) * ps;
2216
2217 }
2218
2219 //==========================================================================
2220  
2221 } // end namespace Pythia8