sistripvvi::VVIObj Class Reference

#include <VVIObj.h>

Public Member Functions
double	fcn (double x) const
void	limits (double &xl, double &xu) const
	density (mode=0) or distribution (mode=1) function More...
	VVIObj (double kappa=0.01, double beta2=1., int mode=0)
	Constructor. More...

Private Attributes
double	a_ [155]
double	b_ [155]
const int	mode_
	returns the limits on the non-zero (mode=0) or normalized region (mode=1) More...
double	omega_
double	t0_
double	t1_
double	t_
double	x0_

Detailed Description

Port of CERNLIB routines vvidis/vviden (G116) to calculate higher quality Vavilov density and distribution functions

Definition at line 25 of file VVIObj.h.

Constructor & Destructor Documentation

VVIObj()

VVIObj::VVIObj	(	double	kappa = 0.01,
		double	beta2 = 1.,
		int	mode = 0
	)

Constructor.

Constructor Set Vavilov parameters kappa and beta2 and define whether to calculate density fcn or distribution fcn

Parameters

kappa	- (input) Vavilov kappa parameter [0.01 (Landau-like) < kappa < 10. (Gaussian-like)]
beta2	- (input) Vavilov beta2 parameter (square of particle speed in v/c units)
mode	- (input) set to 0 to calculate the density function and to 1 to calculate the distribution function

Definition at line 49 of file VVIObj.cc.

References a_, b_, alignmentValidation::c1, funct::cos(), ztail::d, d1, sistripvvi::VVIObjDetails::dzero(), MillePedeFileConverter_cfg::e, JetChargeProducer_cfi::exp, sistripvvi::VVIObjDetails::expint(), sistripvvi::VVIObjDetails::f1(), sistripvvi::VVIObjDetails::f2(), dqmdumpme::k, hgcalLayerClusters_cfi::kappa, cmsLHEtoEOSManager::l, dqm-mbProfile::log, mode_, dqmiodumpmetadata::n, omega_, submitPVResolutionJobs::q, funct::sin(), sistripvvi::VVIObjDetails::sincosint(), t0_, t1_, t_, x, x0_, and testProducerWithPsetDescEmpty_cfi::x1.

                                                     : mode_(mode) {
    const double xp[9] = {9.29, 2.47, .89, .36, .15, .07, .03, .02, 0.0};
    const double xq[7] = {.012, .03, .08, .26, .87, 3.83, 11.0};
    double h_[7];
    double q, u, x, c1, c2, c3, c4, d1, h4, h5, h6, q2, x1, d, ll, ul, xf1, xf2, rv;
    int lp, lq, k, l, n;

    // Make sure that the inputs are reasonable

    if (kappa < 0.01)
      kappa = 0.01;
    if (kappa > 10.)
      kappa = 10.;
    if (beta2 < 0.)
      beta2 = 0.;
    if (beta2 > 1.)
      beta2 = 1.;

    h_[4] = 1. - beta2 * 0.42278433999999998 + 7.6 / kappa;
    h_[5] = beta2;
    h_[6] = 1. - beta2;
    h4 = -7.6 / kappa - (beta2 * .57721566 + 1);
    h5 = log(kappa);
    h6 = 1. / kappa;
    t0_ = (h4 - h_[4] * h5 - (h_[4] + beta2) * (log(h_[4]) + VVIObjDetails::expint(h_[4])) + exp(-h_[4])) / h_[4];

    // Set up limits for the root search

    for (lp = 0; lp < 9; ++lp) {
      if (kappa >= xp[lp])
        break;
    }
    ll = -lp - 1.5;
    for (lq = 0; lq < 7; ++lq) {
      if (kappa <= xq[lq])
        break;
    }
    ul = lq - 6.5;
    //  double (*fp2)(double) = reinterpret_cast<double(*)(double)>(&VVIObj::f2);
    VVIObjDetails::dzero(ll, ul, u, rv, 1.e-5, 1000, std::bind(&VVIObjDetails::f2, std::placeholders::_1, h_));
    q = 1. / u;
    t1_ = h4 * q - h5 - (beta2 * q + 1) * (log((fabs(u))) + VVIObjDetails::expint(u)) + exp(-u) * q;
    t_ = t1_ - t0_;
    omega_ = 6.2831853000000004 / t_;
    h_[0] = kappa * (beta2 * .57721566 + 2.) + 9.9166128600000008;
    if (kappa >= .07) {
      h_[0] += 6.90775527;
    }
    h_[1] = beta2 * kappa;
    h_[2] = h6 * omega_;
    h_[3] = omega_ * 1.5707963250000001;
    //  double (*fp1)(double) = reinterpret_cast<double(*)(double)>(&VVIObj::f1);
    VVIObjDetails::dzero(5., 155., x0_, rv, 1.e-5, 1000, std::bind(&VVIObjDetails::f1, std::placeholders::_1, h_));
    n = x0_ + 1.;
    d = exp(kappa * (beta2 * (.57721566 - h5) + 1.)) * .31830988654751274;
    a_[n - 1] = 0.;
    if (mode_ == 0) {
      a_[n - 1] = omega_ * .31830988654751274;
    }
    q = -1.;
    q2 = 2.;
    for (k = 1; k < n; ++k) {
      l = n - k;
      x = omega_ * k;
      x1 = h6 * x;
      VVIObjDetails::sincosint(x1, c2, c1);
      c1 = log(x) - c1;
      c3 = sin(x1);
      c4 = cos(x1);
      xf1 = kappa * (beta2 * c1 - c4) - x * c2;
      xf2 = x * c1 + kappa * (c3 + beta2 * c2) + t0_ * x;
      if (mode_ == 0) {
        d1 = q * d * omega_ * exp(xf1);
        a_[l - 1] = d1 * cos(xf2);
        b_[l - 1] = -d1 * sin(xf2);
      } else {
        d1 = q * d * exp(xf1) / k;
        a_[l - 1] = d1 * sin(xf2);
        b_[l - 1] = d1 * cos(xf2);
        a_[n - 1] += q2 * a_[l - 1];
      }
      q = -q;
      q2 = -q2;
    }

  }  // VVIObj

Member Function Documentation

fcn()

double VVIObj::fcn

(

double

x

)

const

Vavilov function method Returns density fcn (mode=0) or distribution fcn (mode=1)

Parameters

x	- (input) Argument of function [typically defined as (Q-mpv)/sigma]

Definition at line 142 of file VVIObj.cc.

References a0, testProducerWithPsetDescEmpty_cfi::a2, a_, b0, b1, b2, b_, funct::cos(), f, dqmdumpme::k, mode_, dqmiodumpmetadata::n, omega_, funct::sin(), t0_, t1_, t_, x, and x0_.

                                   {
    // Local variables

    double f, u, y, a0, a1;
    double a2 = 0.;
    double b1, b0, b2, cof;
    int k, n, n1;

    n = x0_;
    if (x < t0_) {
      f = 0.;
    } else if (x <= t1_) {
      y = x - t0_;
      u = omega_ * y - 3.141592653589793;
      cof = cos(u) * 2.;
      a1 = 0.;
      a0 = a_[0];
      n1 = n + 1;
      for (k = 2; k <= n1; ++k) {
        a2 = a1;
        a1 = a0;
        a0 = a_[k - 1] + cof * a1 - a2;
      }
      b1 = 0.;
      b0 = b_[0];
      for (k = 2; k <= n; ++k) {
        b2 = b1;
        b1 = b0;
        b0 = b_[k - 1] + cof * b1 - b2;
      }
      f = (a0 - a2) * .5 + b0 * sin(u);
      if (mode_ != 0) {
        f += y / t_;
      }
    } else {
      f = 0.;
      if (mode_ != 0) {
        f = 1.;
      }
    }
    return f;
  }  // fcn

limits()

void VVIObj::limits	(	double &	xl,
		double &	xu
	)		const

density (mode=0) or distribution (mode=1) function

Vavilov limits method

Parameters

xl	- (output) Smallest value of the argument for the density and the beginning of the normalized region for the distribution
xu	- (output) Largest value of the argument for the density and the end of the normalized region for the distribution

Definition at line 191 of file VVIObj.cc.

References t0_, and t1_.

                                                  {
    xl = t0_;
    xu = t1_;
    return;
  }  // limits

Member Data Documentation

a_

double sistripvvi::VVIObj::a_[155]

private

Definition at line 42 of file VVIObj.h.

Referenced by fcn(), and VVIObj().

b_

double sistripvvi::VVIObj::b_[155]

private

Definition at line 43 of file VVIObj.h.

Referenced by fcn(), and VVIObj().

mode_

const int sistripvvi::VVIObj::mode_

private

returns the limits on the non-zero (mode=0) or normalized region (mode=1)

set to 0 to calculate the density function and to 1 to calculate the distribution function

Definition at line 36 of file VVIObj.h.

Referenced by fcn(), and VVIObj().

omega_

double sistripvvi::VVIObj::omega_

private

Definition at line 40 of file VVIObj.h.

Referenced by fcn(), and VVIObj().

t0_

double sistripvvi::VVIObj::t0_

private

Definition at line 37 of file VVIObj.h.

Referenced by fcn(), limits(), and VVIObj().

t1_

double sistripvvi::VVIObj::t1_

private

Definition at line 38 of file VVIObj.h.

Referenced by fcn(), limits(), and VVIObj().

t_

double sistripvvi::VVIObj::t_

private

Definition at line 39 of file VVIObj.h.

Referenced by fcn(), and VVIObj().

x0_

double sistripvvi::VVIObj::x0_

private

Definition at line 41 of file VVIObj.h.

Referenced by fcn(), and VVIObj().