source: HiSusy/trunk/Pythia8Sgluon/Sigma2SgluonSgluonBase.cxx @ 8

#include "Sigma2SgluonSgluonBase.h"
#include <iostream>
#include <cmath>

Sigma2SgluonSgluonBase::Sigma2SgluonSgluonBase(){
}
Sigma2SgluonSgluonBase::Sigma2SgluonSgluonBase(int idResonanceIn): m_idsgluon(idResonanceIn), m_idsgluonBar(-idResonanceIn){
}
Sigma2SgluonSgluonBase::~Sigma2SgluonSgluonBase(){
}
double Sigma2SgluonSgluonBase::dsigmadcostheta(double ctp){
  ctp=0;
  std::cout << "ERROR YOU SHOULD NOT HAVE LANDED HERE: Sigma2SgluonSgluonBase::dsigmadcostheta " << std::endl;
  return ctp;
}
double Sigma2SgluonSgluonBase::convFactordsigmadcostheta2dsigmadt(double sHat, double mass) {
  return  sHat/2 * sqrtf(1.-4.*mass*mass/sHat);
}
double Sigma2SgluonSgluonBase::integrate(double A, double B, double EPS) {
  double AA,BB,C1,C2,U,S8,S16;
  //  bool MFLAG,RFLAG;
  /*
    C     ******************************************************************
    C
    C     ADAPTIVE DOUBLE PRECISION GAUSSIAN QUADRATURE.
    C
    C     DGAUSS IS SET EQUAL TO THE APPROXIMATE VALUE OF THE INTEGRAL OF
    C     THE FUNCTION F OVER THE INTERVAL (A,B), WITH ACCURACY PARAMETER
    C     EPS.
    C
    C     ******************************************************************
    C
  */
  const double W[12]={ 0.1012285362903762591525313543e0,
                       0.2223810344533744705443559944e0,
                       0.3137066458778872873379622020e0,
                       0.3626837833783619829651504493e0,
                       0.2715245941175409485178057246e-1,
                       0.6225352393864789286284383699e-1,
                       0.9515851168249278480992510760e-1,
                       0.1246289712555338720524762822e0,
                       0.1495959888165767320815017305e0,
                       0.1691565193950025381893120790e0,
                       0.1826034150449235888667636680e0,
                       0.1894506104550684962853967232e0};
  
  const double X[12]={ 0.9602898564975362316835608686e0,
                       0.7966664774136267395915539365e0,
                       0.5255324099163289858177390492e0,
                       0.1834346424956498049394761424e0,
                       0.9894009349916499325961541735e0,
                       0.9445750230732325760779884155e0,
                       0.8656312023878317438804678977e0,
                       0.7554044083550030338951011948e0,
                       0.6178762444026437484466717640e0,
                       0.4580167776572273863424194430e0,
                       0.2816035507792589132304605015e0,
                       0.9501250983763744018531933543e-1};


  //  START.
  double result = 0.0e0;
  if (B == A) return result;
  double CONST = 0.005e0/(B-A);
  BB=A;
  //
  //  COMPUTATIONAL LOOP
 ONE: 
  AA=BB;
  BB=B;

 TWO:    
  C1=0.5e0*(BB+AA);
  C2=0.5e0*(BB-AA);
  S8=0.0e0;

  for (unsigned int i=0; i<4; i++) {
    U=C2*X[i];
    S8=S8+W[i]*(dsigmadcostheta(C1+U)+dsigmadcostheta(C1-U));
  }

  S8=C2*S8;
  S16=0.0e0;

  for (unsigned int i=4; i<12; i++) {
    U=C2*X[i];
    S16=S16+W[i]*(dsigmadcostheta(C1+U)+dsigmadcostheta(C1-U));
  }

  S16=C2*S16;
  if ( fabs(S16-S8) <= EPS*(1.+fabs(S16)) ) goto FIVE;
  BB=C1;
  if ( 1.e0+fabs(CONST*C2) != 1.e0) goto TWO;
  result=0.0e0;

  //  KERMTR("D103.1",LGFILE,MFLAG,RFLAG);
  //  if (MFLAG) {
  std::cout << "FUNCTION DGAUSS ... TOO HIGH ACCURACY REQUIRED" <<std::endl;
    //  }
    //  if ( !RFLAG) abend();
  return result;
 
 FIVE: 
  result=result+S16;
  if (BB != B) goto ONE;
  return result;
}       
//--------------------------------------------------------------------------