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

New function, a punctured smeared exp

parent 772ebc28
// -*- C++ -*-
// $Id: PuncturedSmearedExp.hh,v 1.1 2004/02/05 15:35:25 boudreau Exp $
// ------------------------------------------------------------------------------//
// This function-object makes an exponential with acceptance holes ("punctures") //
// smeared by a resolution function. //
// //
// //
// Joe Boudreau. //
// //
// ------------------------------------------------------------------------------//
#ifndef _PuncturedSmearedExp_h_
#define _PuncturedSmearedExp_h_
#include "CLHEP/GenericFunctions/AbsFunction.hh"
#include "CLHEP/GenericFunctions/Parameter.hh"
namespace Genfun {
/**
* @author
* @ingroup genfun
*/
class PuncturedSmearedExp: public AbsFunction {
FUNCTION_OBJECT_DEF(PuncturedSmearedExp)
public:
// Constructor
PuncturedSmearedExp();
// Copy constructor
PuncturedSmearedExp(const PuncturedSmearedExp &right);
// Destructor:
virtual ~PuncturedSmearedExp();
// Retreive function value
virtual double operator ()(double argument) const;
virtual double operator ()(const Argument & a) const {return operator() (a[0]);}
// Lifetime of exponential:
Parameter & lifetime();
const Parameter & lifetime() const;
// Width of the gaussian:
Parameter & sigma();
const Parameter & sigma() const;
// Puncture this thing:
void puncture(double min, double max);
// Get the puncture parameters:
Parameter & min(unsigned int i);
Parameter & max(unsigned int i);
const Parameter & min(unsigned int i) const;
const Parameter & max(unsigned int i) const;
private:
// These are for calculating mixing terms.
double pow(double x, int n) const ;
double erfc(double x) const ;
// It is illegal to assign an adjustable constant
const PuncturedSmearedExp & operator=(const PuncturedSmearedExp &right);
Parameter _lifetime;
Parameter _sigma;
std::vector<Parameter> _punctures;
};
} // namespace Genfun
#endif
......@@ -63,6 +63,7 @@ libCLHEP_GenericFunctions_@VERSION@_la_SOURCES = \
Power.cc \
Psi2Hydrogen.cc \
PtRelFcn.cc \
PuncturedSmearedExp.cc \
RCBase.cc \
Rectangular.cc \
ReverseExponential.cc \
......
// -*- C++ -*-
// $Id: PuncturedSmearedExp.cc,v 1.1 2004/02/05 15:35:25 boudreau Exp $
#include "CLHEP/GenericFunctions/PuncturedSmearedExp.hh"
#include <sstream>
namespace Genfun {
FUNCTION_OBJECT_IMP(PuncturedSmearedExp)
PuncturedSmearedExp::PuncturedSmearedExp() :
_lifetime ("Lifetime", 1.0, 0.0), // Bounded from below by zero, by default
_sigma ("Sigma", 1.0, 0.0) // Bounded from below by zero, by default
{
}
PuncturedSmearedExp::PuncturedSmearedExp(const PuncturedSmearedExp & right):
_lifetime (right._lifetime),
_sigma (right._sigma),
_punctures (right._punctures)
{
}
void PuncturedSmearedExp::puncture(double min, double max) {
std::ostringstream mn, mx;
mn << "Min_" << _punctures.size()/2;
mx << "Max_" << _punctures.size()/2;
_punctures.push_back(Parameter(mn.str(), min, 0.0, 10.0));
_punctures.push_back(Parameter(mx.str(), max, 0.0, 10.0));
}
Parameter & PuncturedSmearedExp::min(unsigned int i) {
return _punctures[2*i];
}
const Parameter & PuncturedSmearedExp::min(unsigned int i) const {
return _punctures[2*i];
}
Parameter & PuncturedSmearedExp::max(unsigned int i) {
return _punctures[2*i+1];
}
const Parameter & PuncturedSmearedExp::max(unsigned int i) const {
return _punctures[2*i+1];
}
PuncturedSmearedExp::~PuncturedSmearedExp() {
}
Parameter & PuncturedSmearedExp::lifetime() {
return _lifetime;
}
const Parameter & PuncturedSmearedExp::lifetime() const {
return _lifetime;
}
Parameter & PuncturedSmearedExp::sigma() {
return _sigma;
}
const Parameter & PuncturedSmearedExp::sigma() const {
return _sigma;
}
double PuncturedSmearedExp::operator() (double argument) const {
// Fetch the paramters. This operator does not convolve numerically.
static const double sqrtTwo = sqrt(2.0);
double sigma = _sigma.getValue();
double tau = _lifetime.getValue();
double x = argument;
// Copy:
std::vector<double> punctures(_punctures.size());
for (size_t i=0;i<_punctures.size();i++) punctures[i]=_punctures[i].getValue();
// Overlap elimination:
bool overlap=true;
while (overlap) {
overlap=false;
for (size_t i=0;i<punctures.size()/2;i++) {
std::sort(punctures.begin()+2*i, punctures.begin()+2*i+2);
double min1=punctures[2*i];
double max1=punctures[2*i+1];
for (size_t j=i+1;j<punctures.size()/2;j++) {
std::sort(punctures.begin()+2*j, punctures.begin()+2*j+2);
double min2=punctures[2*j];
double max2=punctures[2*j+1];
if ((min2>min1 && max1>min2) || (min1>min2 && max2<min1)) {
punctures[2*i] =std::min(min1,min2);
punctures[2*i+1]=std::max(max1,max2);
std::vector<double>::iterator t0=punctures.begin()+2*j, t1=t0+2;
punctures.erase(t0,t1);
overlap=true;
break;
}
}
if (overlap) break;
}
}
double expG=0,norm=0;
for (size_t i=0;i<punctures.size()/2;i++) {
double a = punctures[2*i];
double b = punctures[2*i+1];
double alpha = (a/sigma + sigma/tau)/sqrtTwo;
double beta = (b/sigma + sigma/tau)/sqrtTwo;
double delta = 1/sqrtTwo/sigma;
norm += 2*tau*exp(1/(4*delta*delta*tau*tau))*(exp(-alpha/(delta*tau)) - exp(-beta/(delta*tau)));
expG += ((erfc(alpha-delta*x) - erfc(beta-delta*x))*exp(-x/tau) );
}
// protection:
if (norm==0) return norm;
return expG/norm;
}
double PuncturedSmearedExp::erfc(double x) const {
// This is accurate to 7 places.
// See Numerical Recipes P 221
double t, z, ans;
z = (x < 0) ? -x : x;
t = 1.0/(1.0+.5*z);
ans = t*exp(-z*z-1.26551223+t*(1.00002368+t*(0.37409196+t*(0.09678418+
t*(-0.18628806+t*(0.27886807+t*(-1.13520398+t*(1.48851587+
t*(-0.82215223+t*0.17087277 ))) ))) )));
if ( x < 0 ) ans = 2.0 - ans;
return ans;
}
double PuncturedSmearedExp::pow(double x,int n) const {
double val=1;
for(int i=0;i<n;i++){
val=val*x;
}
return val;
}
} // 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