VVIObjF Class Reference

#include <VVIObjF.h>

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

Private Attributes
float	a_ [155]
float	b_ [155]
const int	mode_
	returns the limits on the non-zero (mode=0) or normalized region (mode=1) More...
float	omega_
float	t0_
float	t1_
float	t_
float	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 VVIObjF.h.

Constructor & Destructor Documentation

VVIObjF() [1/2]

VVIObjF::VVIObjF	(	float	kappa,
		float	beta2,
		int	mode
	)

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 43 of file VVIObjF.cc.

References a_, funct::abs(), b_, DummyCfis::c, alignmentValidation::c1, ztail::d, d1, VVIObjFDetails::dzero(), MillePedeFileConverter_cfg::e, VVIObjFDetails::expint(), f, validate-o2o-wbm::f1, validate-o2o-wbm::f2, myMath::fast_expf(), myMath::fast_logf(), ALCARECOEcalPhiSym_cff::float, isotrackApplyRegressor::k, hltHgcalLayerClustersEE_cfi::kappa, MainPageGenerator::l, mode_, create_idmaps::n, omega_, submitPVResolutionJobs::q, alignCSCRings::s, VVIObjFDetails::sincosint(), t0_, t1_, t_, x, x0_, and testProducerWithPsetDescEmpty_cfi::x1.

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

  // Make sure that the inputs are reasonable

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

  float invKappa = 1.f / kappa;
  h_[4] = 1.f - beta2 * 0.42278433999999998f + (7.6f * invKappa);
  h_[5] = beta2;
  h_[6] = 1.f - beta2;
  h4 = -(7.6f * invKappa) - (beta2 * .57721566f + 1.f);
  h5 = vdt::fast_logf(kappa);
  h6 = invKappa;
  t0_ = (h4 - h_[4] * h5 - (h_[4] + beta2) * (vdt::fast_logf(h_[4]) + VVIObjFDetails::expint(h_[4])) +
         vdt::fast_expf(-h_[4])) /
        h_[4];

  // Set up limits for the root search

  for (lp = 0; lp < 9; ++lp) {
    if (kappa >= xp[lp])
      break;
  }
  ll = -float(lp) - 1.5f;
  for (lq = 0; lq < 7; ++lq) {
    if (kappa <= xq[lq])
      break;
  }
  ul = lq - 6.5f;
  auto f2 = [h_](float x) {
    return h_[4] - x + h_[5] * (vdt::fast_logf(std::abs(x)) + VVIObjFDetails::expint(x)) - h_[6] * vdt::fast_expf(-x);
  };
  VVIObjFDetails::dzero(ll, ul, u, rv, 1.e-3f, 100, f2);
  q = 1. / u;
  t1_ = h4 * q - h5 - (beta2 * q + 1.f) * (vdt::fast_logf((fabs(u))) + VVIObjFDetails::expint(u)) +
        vdt::fast_expf(-u) * q;
  t_ = t1_ - t0_;
  omega_ = 6.2831853000000004f / t_;
  h_[0] = kappa * (beta2 * .57721566f + 2.f) + 9.9166128600000008f;
  if (kappa >= .07) {
    h_[0] += 6.90775527f;
  }
  h_[1] = beta2 * kappa;
  h_[2] = h6 * omega_;
  h_[3] = omega_ * 1.5707963250000001f;
  auto f1 = [h_](float x) { return h_[0] + h_[1] * vdt::fast_logf(h_[2] * x) - h_[3] * x; };
  VVIObjFDetails::dzero(5.f, 155.f, x0_, rv, 1.e-3f, 100, f1);
  n = x0_ + 1.;
  d = vdt::fast_expf(kappa * (beta2 * (.57721566f - h5) + 1.f)) * .31830988654751274f;
  a_[n - 1] = 0.f;
  if (mode_ == 0) {
    a_[n - 1] = omega_ * .31830988654751274f;
  }
  q = -1.;
  q2 = 2.;

  float x[n];
  x[0] = 0.f;
  float x1[n];
  x1[0] = 0.f;
  float c1[n];
  c1[0] = 0.f;
  float c2[n];
  c2[0] = 0.f;
  float c3[n];
  c3[0] = 0.f;
  float c4[n];
  c4[0] = 0.f;
  float s[n];
  s[0] = 0.f;
  float c[n];
  c[0] = 0.f;
  float xf1[n];
  xf1[0] = 0.f;
  float xf2[n];
  xf2[0] = 0.f;
  for (k = 1; k < n; ++k) {
    x[k] = omega_ * k;
    x1[k] = h6 * x[k];
  }
  for (k = 1; k < n; ++k) {
    VVIObjFDetails::sincosint(x1[k], c2[k], c1[k]);
  }
  for (k = 1; k < n; ++k) {
    c1[k] = vdt::fast_logf(x[k]) - c1[k];
  }
  for (k = 1; k < n; ++k) {
    vdt::fast_sincosf(x1[k], c3[k], c4[k]);
    xf1[k] = kappa * (beta2 * c1[k] - c4[k]) - x[k] * c2[k];
    xf2[k] = x[k] * c1[k] + kappa * (c3[k] + beta2 * c2[k]) + t0_ * x[k];
  }
  for (k = 1; k < n; ++k) {
    vdt::fast_sincosf(xf2[k], s[k], c[k]);
  }
  float d1[n];
  d1[0] = 0.f;
  if (mode_ == 0) {
    for (k = 1; k < n; ++k) {
      d1[k] = q * d * omega_ * vdt::fast_expf(xf1[k]);
      q = -q;
    }
    for (k = 1; k < n; ++k) {
      l = n - k;
      a_[l - 1] = d1[k] * c[k];
      b_[l - 1] = -d1[k] * s[k];
    }
  } else {
    for (k = 1; k < n; ++k) {
      d1[k] = q * d * vdt::fast_expf(xf1[k]) / k;
      q = -q;
    }
    for (k = 1; k < n; ++k) {
      l = n - k;
      a_[l - 1] = d1[k] * s[k];
    }
    for (k = 1; k < n; ++k) {
      l = n - k;
      b_[l - 1] = d1[k] * c[k];
    }
    for (k = 1; k < n; ++k) {
      l = n - k;
      a_[n - 1] += q2 * a_[l - 1];
      q2 = -q2;
    }
  }
}  // VVIObjF

VVIObjF() [2/2]

VVIObjF::VVIObjF

(

float

kappa

)

Constructor with kappa only.

Constructor Set Vavilov parameter kappa and calculate the distribution fcn

Parameters

kappa

- (input) Vavilov kappa parameter [0.01 (Landau-like) < kappa < 10. (Gaussian-like)]

Definition at line 189 of file VVIObjF.cc.

References a_, funct::abs(), b_, DummyCfis::c, alignmentValidation::c1, ztail::d, d1, VVIObjFDetails::dzero(), MillePedeFileConverter_cfg::e, VVIObjFDetails::expint(), f, validate-o2o-wbm::f1, validate-o2o-wbm::f2, myMath::fast_expf(), myMath::fast_logf(), ALCARECOEcalPhiSym_cff::float, isotrackApplyRegressor::k, hltHgcalLayerClustersEE_cfi::kappa, MainPageGenerator::l, create_idmaps::n, omega_, submitPVResolutionJobs::q, alignCSCRings::s, VVIObjFDetails::sincosint(), t0_, t1_, t_, x, x0_, and testProducerWithPsetDescEmpty_cfi::x1.

                            : mode_(1) {
  const float xp[9] = {9.29, 2.47, .89, .36, .15, .07, .03, .02, 0.0};
  const float xq[7] = {.012, .03, .08, .26, .87, 3.83, 11.0};
  float h_[5];
  float q, u, h4, h5, h6, q2, d, ll, ul, rv;
  int lp, lq, k, l, n;

  // Make sure that the inputs are reasonable

  if (kappa < 0.01f)
    kappa = 0.01f;
  if (kappa > 10.f)
    kappa = 10.f;

  float invKappa = 1.f / kappa;
  h_[4] = 0.57721566f + (7.6f * invKappa);
  h4 = -(7.6f * invKappa) - 1.57721566f;
  h5 = vdt::fast_logf(kappa);
  h6 = invKappa;
  t0_ = (h4 - h_[4] * h5 - (h_[4] + 1.f) * (vdt::fast_logf(h_[4]) + VVIObjFDetails::expint(h_[4])) +
         vdt::fast_expf(-h_[4])) /
        h_[4];

  // Set up limits for the root search

  for (lp = 0; lp < 9; ++lp) {
    if (kappa >= xp[lp])
      break;
  }
  ll = -float(lp) - 1.5f;
  for (lq = 0; lq < 7; ++lq) {
    if (kappa <= xq[lq])
      break;
  }
  ul = lq - 6.5f;
  auto f2 = [h_](float x) { return h_[4] - x + (vdt::fast_logf(std::abs(x)) + VVIObjFDetails::expint(x)); };
  VVIObjFDetails::dzero(ll, ul, u, rv, 1.e-3f, 100, f2);
  q = 1. / u;
  t1_ = h4 * q - h5 - (q + 1.f) * (vdt::fast_logf((fabs(u))) + VVIObjFDetails::expint(u)) + vdt::fast_expf(-u) * q;
  t_ = t1_ - t0_;
  omega_ = 6.2831853000000004f / t_;
  h_[0] = kappa * 2.57721566f + 9.9166128600000008f;
  if (kappa >= .07) {
    h_[0] += 6.90775527f;
  }
  h_[1] = kappa;
  h_[2] = h6 * omega_;
  h_[3] = omega_ * 1.5707963250000001f;
  auto f1 = [h_](float x) { return h_[0] + h_[1] * vdt::fast_logf(h_[2] * x) - h_[3] * x; };
  VVIObjFDetails::dzero(5.f, 155.f, x0_, rv, 1.e-3f, 100, f1);
  n = x0_ + 1.;
  d = vdt::fast_expf(kappa * ((0.57721566f - h5) + 1.f)) * .31830988654751274f;
  a_[n - 1] = 0.f;
  q = -1.;
  q2 = 2.;

  float x[n];
  x[0] = 0.f;
  float x1[n];
  x1[0] = 0.f;
  float c1[n];
  c1[0] = 0.f;
  float c2[n];
  c2[0] = 0.f;
  float c3[n];
  c3[0] = 0.f;
  float c4[n];
  c4[0] = 0.f;
  float s[n];
  s[0] = 0.f;
  float c[n];
  c[0] = 0.f;
  float xf1[n];
  xf1[0] = 0.f;
  float xf2[n];
  xf2[0] = 0.f;
  for (k = 1; k < n; ++k) {
    x[k] = omega_ * k;
    x1[k] = h6 * x[k];
  }
  for (k = 1; k < n; ++k) {
    VVIObjFDetails::sincosint(x1[k], c2[k], c1[k]);
  }
  for (k = 1; k < n; ++k) {
    c1[k] = vdt::fast_logf(x[k]) - c1[k];
  }
  for (k = 1; k < n; ++k) {
    vdt::fast_sincosf(x1[k], c3[k], c4[k]);
    xf1[k] = kappa * (c1[k] - c4[k]) - x[k] * c2[k];
    xf2[k] = x[k] * c1[k] + kappa * (c3[k] + c2[k]) + t0_ * x[k];
  }
  for (k = 1; k < n; ++k) {
    vdt::fast_sincosf(xf2[k], s[k], c[k]);
  }
  float d1[n];
  d1[0] = 0.f;
  for (k = 1; k < n; ++k) {
    d1[k] = q * d * vdt::fast_expf(xf1[k]) / k;
    q = -q;
  }
  for (k = 1; k < n; ++k) {
    l = n - k;
    a_[l - 1] = d1[k] * s[k];
  }
  for (k = 1; k < n; ++k) {
    l = n - k;
    b_[l - 1] = d1[k] * c[k];
  }
  for (k = 1; k < n; ++k) {
    l = n - k;
    a_[n - 1] += q2 * a_[l - 1];
    q2 = -q2;
  }

}  // VVIObjF with kappa only

Member Function Documentation

fcn()

float VVIObjF::fcn

(

float

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 311 of file VVIObjF.cc.

References a0, isotrackTrainRegressor::a1, isotrackTrainRegressor::a2, a_, b0, b1, b2, b_, f, isotrackApplyRegressor::k, mode_, create_idmaps::n, omega_, t0_, t1_, t_, x, x0_, and y.

Referenced by SiPixelTemplateReco::PixelTempReco1D(), and SiPixelTemplateReco2D::PixelTempReco2D().

                                {
  // Local variables

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

  n = x0_;
  if (x < t0_) {
    f = 0.f;
  } else if (x <= t1_) {
    y = x - t0_;
    u = omega_ * y - 3.141592653589793f;
    float su, cu;
    vdt::fast_sincosf(u, su, cu);
    cof = cu * 2.f;
    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) * .5f + b0 * su;
    if (mode_ != 0) {
      f += y / t_;
    }
  } else {
    f = 0.f;
    if (mode_ != 0) {
      f = 1.f;
    }
  }
  return f;
}  // fcn

limits()

void VVIObjF::limits	(	float &	xl,
		float &	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 362 of file VVIObjF.cc.

References t0_, and t1_.

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

Member Data Documentation

a_

float VVIObjF::a_[155]

private

Definition at line 42 of file VVIObjF.h.

Referenced by fcn(), and VVIObjF().

b_

float VVIObjF::b_[155]

private

Definition at line 43 of file VVIObjF.h.

Referenced by fcn(), and VVIObjF().

mode_

const int VVIObjF::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 VVIObjF.h.

Referenced by fcn(), and VVIObjF().

omega_

float VVIObjF::omega_

private

Definition at line 40 of file VVIObjF.h.

Referenced by fcn(), and VVIObjF().

t0_

float VVIObjF::t0_

private

Definition at line 37 of file VVIObjF.h.

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

t1_

float VVIObjF::t1_

private

Definition at line 38 of file VVIObjF.h.

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

t_

float VVIObjF::t_

private

Definition at line 39 of file VVIObjF.h.

Referenced by fcn(), and VVIObjF().

x0_

float VVIObjF::x0_

private

Definition at line 41 of file VVIObjF.h.

Referenced by fcn(), and VVIObjF().