Warning, /include/Herwig/PDT/ThreeBodyAllOnCalculator.tcc is written in an unsupported language. File is not indexed.
0001 // -*- C++ -*-
0002 //
0003 // ThreeBodyAllOnCalculator.tcc is a part of Herwig - A multi-purpose Monte Carlo event generator
0004 // Copyright (C) 2002-2019 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 //
0010 // This is the implementation of the non-inlined templated member
0011 // functions of the ThreeBodyAllOnCalculator class.
0012 //
0013 using namespace Herwig;
0014
0015 // shift the variables for the outer integrand and give limits for the inner one
0016 template <class T>
0017 void ThreeBodyAllOnCalculator<T>::outerVariables(double x, Energy2 & low,
0018 Energy2 & upp) const {
0019 // first convert the value of x into the value of souter
0020 if(_mapping[_thechannel]==0) {
0021 _souter = _channelmass[_thechannel]*(_channelmass[_thechannel]+
0022 _channelwidth[_thechannel]*tan(x));
0023 }
0024 else if(_mapping[_thechannel]==1) {
0025 _souter = sqr(_channelmass[_thechannel])*(1.+1./x);
0026 }
0027 else {
0028 _souter = UnitRemoval::E2 * pow(x,1./(_channelpower[_thechannel]+1.));
0029 }
0030 // now the limits of the inner integral
0031 Energy ea(ZERO),eb(ZERO);
0032 Energy rs=sqrt(_souter);
0033 Energy2 eam2(ZERO),ebm2(ZERO);
0034 switch(_channeltype[_thechannel]) {
0035 case 1:
0036 ea = 0.5*(_souter-_m2[1]+_m2[2])/rs;
0037 eam2 = sqr(ea)-_m2[2];
0038 eb = 0.5*(_m2[0]-_souter-_m2[3])/rs;
0039 ebm2 = sqr(eb)-_m2[3];
0040 break;
0041 case 2:
0042 ea = 0.5*(_souter-_m2[1]+_m2[3])/rs;
0043 eam2 = sqr(ea)-_m2[3];
0044 eb = 0.5*(_m2[0]-_souter-_m2[2])/rs;
0045 ebm2 = sqr(eb)-_m2[2];
0046 break;
0047 case 3:
0048 ea = 0.5*(_souter-_m2[2]+_m2[3])/rs;
0049 eam2 = sqr(ea)-_m2[3];
0050 eb = 0.5*(_m2[0]-_souter-_m2[1])/rs;
0051 ebm2 = sqr(eb)-_m2[1];
0052 break;
0053 default:
0054 assert(false);
0055 }
0056 Energy eam = sqrt(max(ZERO,eam2));
0057 Energy ebm = sqrt(max(ZERO,ebm2));
0058 Energy2 sum = sqr(ea+eb);
0059 // calculate the limits
0060 low = sum - sqr(eam+ebm);
0061 upp = sum - sqr(eam-ebm);
0062 }
0063
0064 template <class T>
0065 Energy2 ThreeBodyAllOnCalculator<T>::operator ()(Energy2 y) const {
0066 assert(!std::isnan(double(y/MeV2)));
0067 // set up the values of the s variables
0068 Energy2 s12(ZERO),s23(ZERO),s13(ZERO),
0069 m2sum(_m2[0]+_m2[1]+_m2[2]+_m2[3]);
0070 switch(_channeltype[_thechannel]) {
0071 case 1:
0072 s12 = _souter;
0073 s23 = y;
0074 s13 = m2sum-s12-s23;
0075 break;
0076 case 2:
0077 s23 = y;
0078 s13 = _souter;
0079 s12 = m2sum-s23-s13;
0080 break;
0081 case 3:
0082 s23 = _souter;
0083 s13 = y;
0084 s12 = m2sum-s23-s13;
0085 break;
0086 }
0087 // compute the jacobian
0088 // computer the denominator for the jacobian
0089 InvEnergy2 jacdem = ZERO;
0090 Energy2 sjac(ZERO);
0091 Energy2 rm2,rw2;
0092 for(unsigned int ix=0,N=_channeltype.size(); ix<N; ++ix) {
0093 switch(_channeltype[ix]) {
0094 case 1:
0095 sjac = s12;
0096 break;
0097 case 2:
0098 sjac = s13;
0099 break;
0100 case 3:
0101 sjac = s23;
0102 break;
0103 }
0104 assert(!std::isnan(double(sjac/MeV2)));
0105 InvEnergy2 term;
0106
0107 if(_mapping[ix]==0) {
0108 rm2 = sqr(_channelmass[ix]);
0109 rw2 = sqr(_channelwidth[ix]);
0110 Energy4 tmp = sqr(sjac-rm2) + rw2*rm2;
0111 term = _channelweights[ix]*_channelmass[ix]*_channelwidth[ix]/tmp;
0112 }
0113 else if(_mapping[ix]==1) {
0114 term = _channelweights[ix]*sqr(_channelmass[ix]/(sjac-sqr(_channelmass[ix])));
0115 }
0116 else if(_mapping[ix]==2) {
0117 term = UnitRemoval::InvE2 * _channelweights[ix]*(_channelpower[ix]+1.)*
0118 pow(sjac*UnitRemoval::InvE2, _channelpower[ix]);
0119 }
0120 else
0121 assert(false);
0122 jacdem += term;
0123 }
0124 // now computer the matrix element
0125 return _theME.threeBodyMatrixElement(_mode,_m2[0],s12,s13,
0126 s23,_m[1],_m[2],_m[3])/jacdem;
0127 }
0128
0129 // calculate the width for a given mass
0130 template <class T>
0131 Energy ThreeBodyAllOnCalculator<T>::partialWidth(Energy2 q2) const {
0132 Outer outer(this,_relerr);
0133 _m[0] = sqrt(q2);
0134 _m2[0]=q2;
0135 // check the decay is kinematically allowed
0136 if(_m[0]<_m[1]+_m[2]+_m[3]) return ZERO;
0137 // set up for the different channels
0138 unsigned int N = _channeltype.size();
0139 vector<double> rupp(N,0.),rlow(N,0.);
0140 for(unsigned int ix=0; ix<N; ++ix) {
0141 Energy2 upp(ZERO),low(ZERO);
0142 // work out the kinematic limits
0143 switch(_channeltype[ix]) {
0144 case 1:
0145 upp = sqr(_m[0]-_m[3]);
0146 low = sqr(_m[1]+_m[2]);
0147 break;
0148 case 2:
0149 upp = sqr(_m[0]-_m[2]);
0150 low = sqr(_m[1]+_m[3]);
0151 break;
0152 case 3:
0153 upp = sqr(_m[0]-_m[1]);
0154 low = sqr(_m[2]+_m[3]);
0155 break;
0156 default:
0157 assert(false);
0158 }
0159 // transform them
0160 if(_channelmass[ix] > ZERO) {
0161 if(_channelwidth[ix] > 1e-8*MeV) {
0162 rupp[ix] = atan2((upp-_channelmass[ix]*_channelmass[ix]),
0163 _channelmass[ix]*_channelwidth[ix]);
0164 rlow[ix] = atan2((low-_channelmass[ix]*_channelmass[ix]),
0165 _channelmass[ix]*_channelwidth[ix]);
0166 _mapping[ix] = 0;
0167 if(rupp[ix]/rlow[ix]>0.&&_channelwidth[ix]/_channelmass[ix]<1e-6) {
0168 _mapping[ix] = 1;
0169 Energy2 m2=sqr(_channelmass[ix]);
0170 rupp[ix] = m2/(low-m2);
0171 rlow[ix] = m2/(upp-m2);
0172 }
0173 }
0174 else {
0175 _mapping[ix] = 1;
0176 Energy2 m2=sqr(_channelmass[ix]);
0177 rupp[ix] = m2/(low-m2);
0178 rlow[ix] = m2/(upp-m2);
0179 }
0180 }
0181 else {
0182 _mapping[ix] = 2;
0183 rupp[ix] = pow(upp*UnitRemoval::InvE2, _channelpower[ix]+1.);
0184 rlow[ix] = pow(low*UnitRemoval::InvE2, _channelpower[ix]+1.);
0185 }
0186 }
0187 // perform the integrals for all the different channels
0188 Energy4 sum(ZERO);
0189 for(unsigned int ix=0,N=_channeltype.size(); ix<N; ++ix) {
0190 // perform the integral using GSLIntegrator class
0191 _thechannel=ix;
0192 GSLIntegrator intb(1e-35,_relerr,1000);
0193 sum += _channelweights[ix] * intb.value(outer,rlow[ix],rupp[ix]);
0194 }
0195 // final factors
0196 Energy3 fact = pow<3,1>(Constants::twopi * _m[0]);
0197 return sum/fact/32.;
0198 }
