Warning, /include/Herwig/PDT/ThreeBodyAllOn1IntegralCalculator.tcc is written in an unsupported language. File is not indexed.
0001 // -*- C++ -*-
0002 //
0003 // ThreeBodyAllOn1IntegralCalculator.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 ThreeBodyAllOn1IntegralCalculator class.
0012 //
0013
0014 namespace Herwig{
0015 using namespace ThePEG;
0016
0017 template <class T>
0018 Energy ThreeBodyAllOn1IntegralCalculator<T>::operator() (double x) const {
0019 Energy2 scale;
0020 if(_intmass>ZERO) scale = _intmass*(_intmass+_intwidth*tan(x));
0021 else scale = UnitRemoval::E2 * pow(x,1./(_intpower + 1.));
0022 InvEnergy output=_theDgamma.threeBodydGammads(_mode,_m2[0],scale,_m[1],_m[2],_m[3]);
0023 // the jacobian
0024 InvEnergy2 term;
0025 Energy2 rm2,rw2;
0026 if(_intmass>ZERO) {
0027 rm2 = sqr(_intmass);
0028 rw2 = sqr(_intwidth);
0029 term = _intmass*_intwidth / (sqr(scale-rm2) + rw2 * rm2);
0030 }
0031 else {
0032 term = UnitRemoval::InvE2 * (_intpower+1.)*pow(scale * UnitRemoval::InvE2, _intpower);
0033 }
0034 return output/term;
0035 }
0036
0037 template <class T>
0038 Energy ThreeBodyAllOn1IntegralCalculator<T>::partialWidth(Energy2 scale) const {
0039 _m2[0]=scale;
0040 _m[0]=sqrt(scale);
0041 // limits for the outer integral
0042 Energy2 upp=ZERO,low=ZERO;
0043 switch(_variabletype) {
0044 case 1:
0045 upp = (_m[0]-_m[3])*(_m[0]-_m[3]);
0046 low = (_m[1]+_m[2])*(_m[1]+_m[2]);
0047 break;
0048 case 2:
0049 upp = (_m[0]-_m[2])*(_m[0]-_m[2]);
0050 low = (_m[1]+_m[3])*(_m[1]+_m[3]);
0051 break;
0052 case 3:
0053 upp = (_m[0]-_m[1])*(_m[0]-_m[1]);
0054 low = (_m[2]+_m[3])*(_m[2]+_m[3]);
0055 break;
0056 }
0057 double rupp, rlow;
0058 // transform them
0059 if(_intmass>ZERO) {
0060 rupp = atan2((upp-_intmass*_intmass), _intmass*_intwidth);
0061 rlow = atan2((low-_intmass*_intmass), _intmass*_intwidth);
0062 }
0063 else {
0064 rupp = pow(upp * UnitRemoval::InvE2, _intpower+1.);
0065 rlow = pow(low * UnitRemoval::InvE2, _intpower+1.);
0066 }
0067 return _integrator.value(*this,rlow,rupp);
0068 }
0069
0070 }
0071
