Commit f0cc6c0c authored by Joseph Boudreau's avatar Joseph Boudreau
Browse files

Added three BW functions

parent 74492123
......@@ -75,6 +75,7 @@ pkginclude_HEADERS = \
LogGamma.hh \
LogisticFunction.hh \
Mod.hh \
NonrelativisticBW.hh \
ParameterComposition.hh \
ParameterDifference.hh \
Parameter.hh \
......@@ -91,6 +92,7 @@ pkginclude_HEADERS = \
PuncturedSmearedExp.hh \
RCBase.hh \
Rectangular.hh \
RelativisticBW.hh \
ReverseExponential.hh \
RKIntegrator.hh \
Sigma.hh \
......@@ -106,6 +108,7 @@ pkginclude_HEADERS = \
Tan.hh \
TrivariateGaussian.hh \
Variable.hh \
VoigtProfile.hh \
X.hh \
defs.h
......
// -*- C++ -*-
// $Id:
//---------------------NonrelativisticBWDistribution------------------------//
// //
// //
// Joe Boudreau, June 2011 //
// //
//--------------------------------------------------------------------------//
#ifndef NonrelativisticBWDistribution_h
#define NonrelativisticBWDistribution_h 1
#include "CLHEP/GenericFunctions/AbsFunction.hh"
#include "CLHEP/GenericFunctions/Parameter.hh"
#include "CLHEP/GenericFunctions/IncompleteGamma.hh"
namespace Genfun {
/**
* @author
* @ingroup genfun
*/
class NonrelativisticBWDistribution : public AbsFunction {
FUNCTION_OBJECT_DEF(NonrelativisticBWDistribution)
public:
// Constructor
NonrelativisticBWDistribution();
// Copy constructor
NonrelativisticBWDistribution(const NonrelativisticBWDistribution &right);
// Destructor
virtual ~NonrelativisticBWDistribution();
// Retreive function value
virtual double operator ()(double argument) const;
virtual double operator ()(const Argument & a) const {return operator() (a[0]);}
// Get the paramter alpha
Parameter & mass();
// Get the parameter beta
Parameter & width();
private:
// It is illegal to assign an adjustable constant
const NonrelativisticBWDistribution & operator=(const NonrelativisticBWDistribution &right);
// Here are the two parameters alpha and beta:
Parameter _mass;
Parameter _width;
};
} // namespace Genfun
#endif
// -*- C++ -*-
// $Id:
//---------------------RelativisticBWDistribution---------------------------//
// //
// //
// Joe Boudreau, June 2011 //
// //
//--------------------------------------------------------------------------//
#ifndef RelativisticBWDistribution_h
#define RelativisticBWDistribution_h 1
#include "CLHEP/GenericFunctions/AbsFunction.hh"
#include "CLHEP/GenericFunctions/Parameter.hh"
#include "CLHEP/GenericFunctions/IncompleteGamma.hh"
namespace Genfun {
/**
* @author
* @ingroup genfun
*/
class RelativisticBWDistribution : public AbsFunction {
FUNCTION_OBJECT_DEF(RelativisticBWDistribution)
public:
// Constructor
RelativisticBWDistribution();
// Copy constructor
RelativisticBWDistribution(const RelativisticBWDistribution &right);
// Destructor
virtual ~RelativisticBWDistribution();
// Retreive function value
virtual double operator ()(double argument) const;
virtual double operator ()(const Argument & a) const {return operator() (a[0]);}
// Get the paramter alpha
Parameter & mass();
// Get the parameter beta
Parameter & width();
private:
// It is illegal to assign an adjustable constant
const RelativisticBWDistribution & operator=(const RelativisticBWDistribution &right);
// Here are the two parameters alpha and beta:
Parameter _mass;
Parameter _width;
};
} // namespace Genfun
#endif
// -*- C++ -*-
// $Id:
//---------------------VoigtProfile----------------------------------------//
// //
// //
// Joe Boudreau, June 2011 //
// //
//--------------------------------------------------------------------------//
#ifndef VoigtProfile_h
#define VoigtProfile_h 1
#include "CLHEP/GenericFunctions/AbsFunction.hh"
#include "CLHEP/GenericFunctions/Parameter.hh"
#include "CLHEP/GenericFunctions/IncompleteGamma.hh"
namespace Genfun {
/**
* @author
* @ingroup genfun
*/
class VoigtProfile : public AbsFunction {
FUNCTION_OBJECT_DEF(VoigtProfile)
public:
// Constructor
VoigtProfile();
// Copy constructor
VoigtProfile(const VoigtProfile &right);
// Destructor
virtual ~VoigtProfile();
// Retreive function value
virtual double operator ()(double argument) const;
virtual double operator ()(const Argument & a) const {return operator() (a[0]);}
// Get the paramter alpha
Parameter & mass();
// Get the parameter beta
Parameter & width();
// Get the parameter beta
Parameter & sigma();
private:
// It is illegal to assign an adjustable constant
const VoigtProfile & operator=(const VoigtProfile &right);
// Here are the two parameters alpha and beta:
Parameter _mass;
Parameter _width;
Parameter _sigma;
};
} // namespace Genfun
#endif
......@@ -58,6 +58,7 @@ libCLHEP_GenericFunctions_@VERSION@_a_SOURCES = \
LogGamma.cc \
LogisticFunction.cc \
Mod.cc \
NonrelativisticBW.cc \
Parameter.cc \
ParameterComposition.cc \
ParameterDifference.cc \
......@@ -72,6 +73,7 @@ libCLHEP_GenericFunctions_@VERSION@_a_SOURCES = \
PuncturedSmearedExp.cc \
RCBase.cc \
Rectangular.cc \
RelativisticBW.cc \
ReverseExponential.cc \
RKIntegrator.cc \
Sigma.cc \
......@@ -81,6 +83,7 @@ libCLHEP_GenericFunctions_@VERSION@_a_SOURCES = \
Tan.cc \
TrivariateGaussian.cc \
Variable.cc \
VoigtProfile.cc \
X.cc
libCLHEP_GenericFunctions_@VERSION@_so_OBJECTS = $(patsubst %.cc,$(shareddir)/%.$(OBJEXT),$(libCLHEP_GenericFunctions_@VERSION@_a_SOURCES))
......
#include "CLHEP/GenericFunctions/NonrelativisticBW.hh"
#include "CLHEP/GenericFunctions/Variable.hh"
#include <assert.h>
#include <cmath>
using namespace std;
namespace Genfun {
FUNCTION_OBJECT_IMP(NonrelativisticBWDistribution)
NonrelativisticBWDistribution::NonrelativisticBWDistribution():
_mass("mass", 50, 10, 90),
_width ("width", 5, 0, 100)
{}
NonrelativisticBWDistribution::NonrelativisticBWDistribution(const NonrelativisticBWDistribution & right):
AbsFunction(),
_mass(right._mass),
_width (right._width)
{
}
NonrelativisticBWDistribution::~NonrelativisticBWDistribution() {
}
double NonrelativisticBWDistribution::operator() (double x) const {
double M=_mass.getValue();
double G=_width.getValue()/2.0;
double f = (1.0/M_PI)*G/((x-M)*(x-M) +G*G);
return f;
}
Parameter & NonrelativisticBWDistribution::mass() {
return _mass;
}
Parameter & NonrelativisticBWDistribution::width() {
return _width;
}
} // namespace Genfun
#include "CLHEP/GenericFunctions/RelativisticBW.hh"
#include "CLHEP/GenericFunctions/Variable.hh"
#include <assert.h>
#include <cmath>
using namespace std;
namespace Genfun {
FUNCTION_OBJECT_IMP(RelativisticBWDistribution)
RelativisticBWDistribution::RelativisticBWDistribution():
_mass("mass", 50, 10, 90),
_width ("width", 5, 0, 100)
{}
RelativisticBWDistribution::RelativisticBWDistribution(const RelativisticBWDistribution & right):
AbsFunction(),
_mass(right._mass),
_width (right._width)
{
}
RelativisticBWDistribution::~RelativisticBWDistribution() {
}
double RelativisticBWDistribution::operator() (double x) const {
double M=_mass.getValue();
double G=_width.getValue();
double g=sqrt(M*M*(M*M+G*G));
double k = 2.0*sqrt(2.0)*M*G*g/M_PI/sqrt(M*M+g);
double f = k/((x-M)*(x-M)*(x+M)*(x+M)+M*M*G*G);
return f;
}
Parameter & RelativisticBWDistribution::mass() {
return _mass;
}
Parameter & RelativisticBWDistribution::width() {
return _width;
}
} // namespace Genfun
#include "CLHEP/GenericFunctions/VoigtProfile.hh"
#include "CLHEP/GenericFunctions/Variable.hh"
#include <assert.h>
#include <cmath>
#include <complex>
using namespace std;
inline double Pow(double x,int n) {
double val=1;
for(int i=0;i<n;i++){
val=val*x;
}
return val;
}
inline std::complex<double> nwwerf(std::complex<double> z) {
std::complex<double> zh,r[38],s,t,v;
const double z1 = 1;
const double hf = z1/2;
const double z10 = 10;
const double c1 = 74/z10;
const double c2 = 83/z10;
const double c3 = z10/32;
const double c4 = 16/z10;
const double c = 1.12837916709551257;
const double p = Pow(2.0*c4,33);
double x=z.real();
double y=z.imag();
double xa=(x >= 0) ? x : -x;
double ya=(y >= 0) ? y : -y;
if(ya < c1 && xa < c2){
zh = std::complex<double>(ya+c4,xa);
r[37]= std::complex<double>(0,0);
// do 1 n = 36,1,-1
for(int n = 36; n>0; n--){
t=zh+double(n)*std::conj(r[n+1]);
r[n]=hf*t/std::norm(t);
}
double xl=p;
s=std::complex<double>(0,0);
// do 2 n = 33,1,-1
for(int k=33; k>0; k--){
xl=c3*xl;
s=r[k]*(s+xl);
}
v=c*s;
}
else{
zh=std::complex<double>(ya,xa);
r[1]=std::complex<double>(0,0);
// do 3 n = 9,1,-1
for(int n=9;n>0;n--){
t=zh+double(n)*std::conj(r[1]);
r[1]=hf*t/std::norm(t);
}
v=c*r[1];
}
if(ya == 0) v= std::complex<double>(exp(-xa*xa),v.imag());
if(y < 0) {
v=2.0*std::exp(std::complex<double>(-xa,-ya)*std::complex<double>(xa,ya))-v;
if(x > 0) v=std::conj(v);
}
else{
if(x < 0) v=std::conj(v);
}
return v;
}
namespace Genfun {
FUNCTION_OBJECT_IMP(VoigtProfile)
VoigtProfile::VoigtProfile():
_mass("mass", 50, 10, 90),
_width ("width", 5, 0, 100),
_sigma("sigma", 5, 0, 100)
{}
VoigtProfile::VoigtProfile(const VoigtProfile & right):
AbsFunction(),
_mass(right._mass),
_width (right._width),
_sigma(right._sigma)
{
}
VoigtProfile::~VoigtProfile() {
}
double VoigtProfile::operator() (double x) const {
double M=_mass.getValue();
double G=_width.getValue()/2.0;
double s=_sigma.getValue();
static const double sqrt2=sqrt(2.0);
static const double sqrt2PI=sqrt(2.0*M_PI);
static const std::complex<double> I(0,1);
std::complex<double> z = ((x-M) + I*G)/sqrt2/s;
double f=nwwerf(z).real()/s/sqrt2PI;
return f;
}
Parameter & VoigtProfile::mass() {
return _mass;
}
Parameter & VoigtProfile::width() {
return _width;
}
Parameter & VoigtProfile::sigma() {
return _sigma;
}
} // namespace Genfun
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment