Herwig/Utilities/GaussianIntegrator.tcc

0001 // -*- C++ -*-
0002 //
0003 // GaussianIntegrator.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 GaussianIntegrator class.
0012 //
0013
0014 namespace Herwig {
0015 using namespace ThePEG;
0016
0017 template <class T>
0018 inline GaussianIntegrator::ValT<T>
0019 GaussianIntegrator::value(const T & function,
0020                           const typename T::ArgType lower,
0021                           const typename T::ArgType upper) const {
0022   typedef typename T::ValType ValType;
0023   typedef typename T::ArgType ArgType;
0024   const ValType ValUnit = TypeTraits<ValType>::baseunit();
0025   const ArgType ArgUnit = TypeTraits<ArgType>::baseunit();
0026
0027   // vector for the limits of the bin
0028   vector<double> lowerlim,upperlim;
0029   // start with the whole interval as 1 bin
0030   lowerlim.push_back(lower/ArgUnit);upperlim.push_back(upper/ArgUnit);
0031   // set the minimum bin width
0032   double xmin=_binwidth*abs(upper-lower)/ArgUnit;
0033   // counters for the number of function evals
0034   int neval=0;
0035   // and number of bad intervals
0036   int nbad=0;
0037   // the output value
0038   double output=0.;
0039   // the loop for the evaluation
0040   double mid,wid; unsigned int ibin,ix=0,iorder;
0041   double testvalue,value,tolerance;
0042   do {
0043     // the bin we are doing (always the last one in the list)
0044     ibin = lowerlim.size()-1;
0045     // midpoint and width of the bin
0046     mid=0.5*(upperlim[ibin]+lowerlim[ibin]);
0047     wid=0.5*(upperlim[ibin]-lowerlim[ibin]);
0048     value=0.;
0049     iorder=0;
0050     // compute a trail value using sixth order GQ
0051     for(ix=0;ix<_weights[0].size();++ix) {
0052       value+=_weights[0][ix]
0053         *( function((mid+wid*_abscissae[0][ix])*ArgUnit)
0054           +function((mid-wid*_abscissae[0][ix])*ArgUnit)
0055            )/ValUnit;
0056       ++neval;
0057       if(neval>_maxeval)
0058         CurrentGenerator::log() << "Error in Gaussian Integrator: Setting to zero"
0059                                 << endl;
0060     }
0061     value *=wid;
0062     // compute more accurate answers using higher order GQ
0063     do {
0064       // use the next order of quadrature
0065       testvalue=value;
0066       ++iorder;
0067       value=0.;
0068       for(ix=0;ix<_weights[iorder].size();++ix) {
0069         value+=_weights[iorder][ix]*
0070           ( function((mid+wid*_abscissae[iorder][ix])*ArgUnit)
0071            +function((mid-wid*_abscissae[iorder][ix])*ArgUnit)
0072             )/ValUnit;
0073         ++neval;
0074         if(neval>_maxeval)
0075            CurrentGenerator::log() << "Error in Gaussian Integrator: Setting to zero"
0076                                    << endl;
0077       }
0078       value *=wid;
0079       tolerance=max(_abserr,_relerr*abs(value));
0080     }
0081     // keep going if possible and not accurate enough
0082     while(iorder<_weights.size()-1&&abs(testvalue-value)>tolerance);
0083     // now decide what to do
0084     // accept this value
0085     if(abs(testvalue-value)<tolerance) {
0086       output+=value;
0087       lowerlim.pop_back();upperlim.pop_back();
0088     }
0089     // bin too small to redivide contribution set to zero
0090     else if(wid<xmin) {
0091       ++nbad;
0092       lowerlim.pop_back(); upperlim.pop_back();
0093     }
0094     // otherwise split the bin into two
0095     else {
0096       // reset the limits for the bin
0097       upperlim[ibin]=mid;
0098       // set up a new bin
0099       lowerlim.push_back(mid);
0100       upperlim.push_back(mid+wid);
0101     }
0102   }
0103   // keep going if there's still some bins to evaluate
0104   while(lowerlim.size()>0);
0105   // output an error message if needed
0106   if(nbad!=0)
0107     CurrentGenerator::log() << "Error in GaussianIntegrator: Bad Convergence for "
0108                             << nbad << "intervals" << endl;
0109   // return the answer
0110   return output * ValUnit * ArgUnit;
0111 }
0112
0113 }