Herwig/Decay/ResonanceHelpers.h

0001 // -*- C++ -*-
0002 //
0003 // ResonanceHelpers.h is a part of Herwig - A multi-purpose Monte Carlo event generator
0004 // Copyright (C) 2018 The Herwig Collaboration
0005 //
0006 // Herwig is licenced under version 3 of the GPL, see COPYING for details.
0007 // Please respect the MCnet academic guidelines, see GUIDELINES for details.
0008 //
0009 #ifndef HERWIG_ResonanceHelpers_H
0010 #define HERWIG_ResonanceHelpers_H
0011
0012 namespace Herwig {
0013 using namespace ThePEG;
0014 namespace Resonance {
0015
0016 /**
0017  *   The velocity squared
0018  */
0019 inline double beta2(const Energy2 & s, const Energy & m1, const Energy & m2) {
0020   return max(0.,(1.-sqr(m1+m2)/s)*(1.-sqr(m1-m2)/s));
0021 }
0022
0023 /**
0024  *   The velocity
0025  */
0026 inline double beta(const Energy2 & s, const Energy & m1, const Energy & m2) {
0027   return sqrt(beta2(s,m1,m2));
0028 }
0029
0030 /**
0031  *  The derivative of the function \f$\hat{H}(s)\f$ for the GS Breit-Wigner evaluated
0032  *  at the resonance mass
0033  */
0034 inline double dHhatds(const Energy & mRes, const Energy & gamma,
0035               const Energy & m1, const Energy & m2) {
0036   double v2 = beta2(sqr(mRes),m1,m2);
0037   double v = sqrt(v2);
0038   double r = (sqr(m1) + sqr(m2))/sqr(mRes);
0039   return gamma/Constants::pi/mRes/v2*
0040     ((3.-2.*v2- 3.*r)*log((1.+v)/(1.-v)) + 2.*v*(1.- r/(1.-v2)));
0041 }
0042
0043 /**
0044  *  The  \f$\hat{H}(s)\f$ function for the GS Breit-Wigner
0045  */
0046 inline Energy2 Hhat(const Energy2 & s, const Energy & mRes, const Energy & gamma,
0047             const Energy & m1, const Energy & m2) {
0048   double vR = beta(sqr(mRes),m1,m2);
0049   double v  = beta(    s    ,m1,m2);
0050   return gamma/mRes/Constants::pi*s*pow(v/vR,3)*log((1.+v)/(1.-v));
0051 }
0052
0053 /**
0054  *  The \f$H(s)\f$ function for the GS Breit-Wigner (input the derivative term)
0055  */
0056 inline Energy2 H(const Energy2 & s, const Energy & mRes, const Energy & gamma,
0057          const Energy & m1, const Energy & m2,
0058                  const double & dH, const Energy2 & Hres) {
0059   if(s!=ZERO)
0060     return Hhat(s,mRes,gamma,m1,m2) - Hres - (s-sqr(mRes))*dH;
0061   else
0062     return -2.*sqr(m1+m2)/Constants::pi*gamma/mRes/pow(beta(sqr(mRes),m1,m2),3) - Hres + sqr(mRes)*dH;
0063 }
0064
0065 /**
0066  *  The \f$H(s)\f$ function for the GS Breit-Wigner (with calc of derivative term)
0067  */
0068 inline Energy2 H(const Energy2 & s, const Energy & mRes, const Energy & gamma,
0069          const Energy & m1, const Energy & m2) {
0070   double dH = dHhatds(mRes,gamma,m1,m2);
0071   Energy2 Hres = Hhat(sqr(mRes),mRes,gamma,m1,m2);
0072   return H(s,mRes,gamma,m1,m2,dH,Hres);
0073 }
0074
0075 /**
0076  *    The \f$p\f$-wave runningwidth
0077  */
0078 inline Energy gammaP(const Energy2 & s, const Energy & mRes, const Energy & gamma,
0079              const Energy & m1, const Energy & m2) {
0080   double v2 = beta2(s,m1,m2);
0081   if(v2<=0.) return ZERO;
0082   double vR2 = beta2(sqr(mRes),m1,m2);
0083   double rp = sqrt(v2/vR2);
0084   return sqrt(s)/mRes*pow(rp,3)*gamma;
0085 }
0086
0087 /**
0088  *    The \f$p\f$-wave runningwidth
0089  */
0090 inline Energy gammaD(const Energy2 & s, const Energy & mRes, const Energy & gamma,
0091              const Energy & m1, const Energy & m2) {
0092   double v2 = beta2(s,m1,m2);
0093   if(v2<=0.) return ZERO;
0094   double vR2 = beta2(sqr(mRes),m1,m2);
0095   double rp = sqrt(v2/vR2);
0096   return pow(sqrt(s)/mRes,3)*pow(rp,5)*gamma;
0097 }
0098
0099 /**
0100  *    The \f$p\f$-wave runningwidth
0101  */
0102 inline Energy gammaS(const Energy2 & s, const Energy & mRes, const Energy & gamma,
0103              const Energy & m1, const Energy & m2) {
0104   double v2 = beta2(s,m1,m2);
0105   if(v2<=0.) return ZERO;
0106   double vR2 = beta2(sqr(mRes),m1,m2);
0107   double rp = sqrt(v2/vR2);
0108   return mRes/sqrt(s)*rp*gamma;
0109 }
0110
0111 /**
0112  *  The GS form of the Breit-Wigner distribution
0113  */
0114 inline Complex BreitWignerGS(const Energy2 & s, const Energy & mRes, const Energy & gamma,
0115                  const Energy & m1, const Energy & m2,
0116                              const Energy2 & H0, const double &dH, const Energy2 & Hres) {
0117   Energy2 mR2=sqr(mRes);
0118   return (mR2+H0)/(mR2-s+H(s,mRes,gamma,m1,m2,dH,Hres)-Complex(0.,1.)*sqrt(s)*gammaP(s,mRes,gamma,m1,m2));
0119 }
0120
0121 /**
0122  *  The GS form of the Breit-Wigner distribution
0123  */
0124 inline Complex BreitWignerGS(const Energy2 & s, const Energy & mRes,
0125                              const Energy & gamma,
0126                  const Energy & m1, const Energy & m2) {
0127   double dH = dHhatds(mRes,gamma,m1,m2);
0128   Energy2 Hres = Hhat(sqr(mRes),mRes,gamma,m1,m2);
0129   Energy2 H0 = H(ZERO,mRes,gamma,m1,m2,dH,Hres);
0130   return BreitWignerGS(s,mRes,gamma,m1,m2,H0,dH,Hres);
0131 }
0132
0133 /**
0134  *  Standard \f$p\f$-wave Breit-Wigner
0135  */
0136 inline Complex BreitWignerPWave(const Energy2 & s, const Energy & mRes, const Energy & gamma,
0137                 const Energy & m1, const Energy & m2) {
0138   Energy2 mR2=sqr(mRes);
0139   return mR2/(mR2-s-Complex(0.,1.)*sqrt(s)*gammaP(s,mRes,gamma,m1,m2));
0140 }
0141
0142 /**
0143  *  Standard \f$s\f$-wave Breit-Wigner
0144  */
0145 inline Complex BreitWignerSWave(const Energy2 & s, const Energy & mRes, const Energy & gamma,
0146                 const Energy & m1, const Energy & m2) {
0147   Energy2 mR2=sqr(mRes);
0148   return mR2/(mR2-s-Complex(0.,1.)*sqrt(s)*gammaS(s,mRes,gamma,m1,m2));
0149 }
0150
0151 /**
0152  *  Standard \f$d\f$-wave Breit-Wigner
0153  */
0154 inline Complex BreitWignerDWave(const Energy2 & s, const Energy & mRes, const Energy & gamma,
0155                 const Energy & m1, const Energy & m2) {
0156   Energy2 mR2=sqr(mRes);
0157   return mR2/(mR2-s-Complex(0.,1.)*sqrt(s)*gammaD(s,mRes,gamma,m1,m2));
0158 }
0159
0160 /**
0161  *  Standard fixed width Breit-Wigner, no width in numerator
0162  */
0163 inline Complex BreitWignerFW(const Energy2 & s, const Energy & mRes, const Energy & gamma) {
0164   Energy2 mR2=sqr(mRes);
0165   return mR2/(mR2-s-Complex(0.,1.)*mRes*gamma);
0166 }
0167
0168 /**
0169  *  Standard fixed width Breit-Wigner, width in numerator
0170  */
0171 inline Complex BreitWignerFW_GN(const Energy2 & s, const Energy & mRes, const Energy & gamma) {
0172   Energy2 mR2=sqr(mRes);
0173   complex<Energy2> fact = mR2 - Complex(0.,1.)*mRes*gamma;
0174   return fact/(fact-s);
0175 }
0176
0177 /**
0178  *   The \f$H\f$ function from 0512180
0179  */
0180 Complex H(const Energy & mass, const Energy & width, const Energy2 & sp, const Energy2 & sm,
0181       const Energy2 & s0, const Energy & mp, const Energy & m0) {
0182   return
0183     Resonance::BreitWignerPWave(sp,mass,width,mp,m0)+
0184     Resonance::BreitWignerPWave(sm,mass,width,mp,m0)+
0185     Resonance::BreitWignerPWave(s0,mass,width,mp,mp);
0186 }
0187
0188 /**
0189  *  Sum over \f$p\f$-wave resonaces
0190  */
0191 template<typename Value>
0192 Complex F_rho(const Energy2 & s,
0193           const vector<Value> weights,
0194           const vector<Energy> & mass,
0195           const vector<Energy> & width,
0196           const Energy & m1, const Energy & m2) {
0197   Value norm(0.);
0198   Complex output;
0199   for(unsigned int ix=0;ix<weights.size();++ix) {
0200     norm += weights[ix];
0201     output += weights[ix]*
0202       BreitWignerPWave(s,mass[ix],width[ix],m1,m2);
0203   }
0204   return output/norm;
0205 }
0206
0207 double ga1(const Energy2 &s) {
0208   static const Energy mpi=0.13957*GeV;
0209   if(s<9.*sqr(mpi)) {
0210     return 0.;
0211   }
0212   else if(s>0.838968432668*GeV2) {
0213     double Q2 = s/GeV2;
0214     return 1.623*Q2+10.38-9.32/Q2+0.65/sqr(Q2);
0215   }
0216   else {
0217     double Q2 = (s-9.*sqr(mpi))/GeV2;
0218     return 4.1*Q2*sqr(Q2)*(1.-3.3*Q2+5.8*sqr(Q2));
0219   }
0220 }
0221
0222 /**
0223  *  GS form of the \f$a_1\f$ Breit-Wigner
0224  */
0225 Complex BreitWignera1(const Energy2 & s, const Energy & mRes,
0226               const Energy & gamma) {
0227   Energy2 mR2 = sqr(mRes);
0228   return mR2/(mR2-s-Complex(0.,1)*gamma*mRes*ga1(s)/ga1(mR2));
0229 }
0230
0231 /**
0232  *  Difference between two resonances
0233  */
0234 complex<InvEnergy2> BreitWignerDiff(const Energy2 & s,
0235                     const Energy & mRes1, const Energy & gamma1,
0236                     const Energy & mRes2, const Energy & gamma2,
0237                     const Energy & m1, const Energy & m2) {
0238   return
0239     BreitWignerPWave(s,mRes1,gamma1,m1,m2)/sqr(mRes1)-
0240     BreitWignerPWave(s,mRes2,gamma2,m1,m2)/sqr(mRes2);
0241 }
0242 }
0243
0244 }
0245 #endif