Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Simplify VVIObjF object so it takes kappa only in the constructor #37533

Merged
merged 2 commits into from
Apr 20, 2022
Merged

Simplify VVIObjF object so it takes kappa only in the constructor #37533

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
dda1383
Simplify VVIObjF object
tvami Apr 11, 2022
4fb03c3
Add a warning to change both VVIObjF constructors
tvami Apr 14, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
5 changes: 3 additions & 2 deletions RecoLocalTracker/SiPixelRecHits/interface/VVIObjF.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,7 @@
// V1.1 - make dzero call both fcns with a switch
// V1.2 - remove inappriate initializers and add methods to return non-zero/normalized region
// V2.0 - restructuring and speed improvements by V. Innocente
// V3.0 - further simplification and speedup by Tamas Vami
//

#ifndef VVIObjF_h
Expand All @@ -23,15 +24,15 @@
// ***********************************************************************************************************************
class VVIObjF {
public:
VVIObjF(float kappa = 0.01, float beta2 = 1., int mode = 0); //!< Constructor
VVIObjF(float kappa, float beta2, int mode); //!< Constructor
VVIObjF(float kappa); //!< Constructor with kappa only

float fcn(float x) const; //! density (mode=0) or distribution (mode=1) function
void limits(float& xl,
float& xu) const; //! returns the limits on the non-zero (mode=0) or normalized region (mode=1)

private:
// Vavilov distribution parameters (inputs and common block /G116C1/)

const int mode_; //!< set to 0 to calculate the density function and to 1 to calculate the distribution function
float t0_;
float t1_;
Expand Down
5 changes: 3 additions & 2 deletions RecoLocalTracker/SiPixelRecHits/src/SiPixelTemplateReco.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -44,10 +44,11 @@
// V9.00 - Set QProb = Q/Q_avg when calcultion is turned off, use fbin definitions of Qbin
// V10.00 - Use new template object to reco Phase 1 FPix hits
// V10.01 - Fix memory overwriting bug
// V10.10 - Change VVIObjF so it only reads kappa
//
//
// Created by Morris Swartz on 10/27/06.
//
// VVIObjF object simplification by Tamas Vami
//

#ifndef SI_PIXEL_TEMPLATE_STANDALONE
Expand Down Expand Up @@ -1181,7 +1182,7 @@ int SiPixelTemplateReco::PixelTempReco1D(int id,
beta2 = 1.;
if (use_VVIObj) {
// VVIObj is a private port of CERNLIB VVIDIS
VVIObjF vvidist(kappa, beta2, 1);
VVIObjF vvidist(kappa);
prvav = vvidist.fcn(xvav);
} else {
// Use faster but less accurate TMath Vavilov distribution function
Expand Down
6 changes: 3 additions & 3 deletions RecoLocalTracker/SiPixelRecHits/src/SiPixelTemplateReco2D.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -14,11 +14,12 @@
// 2.70 - Change convergence criterion to require it in both planes [it was either]
// 2.80 - Change 3D to 2D
// 2.90 - Fix divide by zero for separate 1D convergence branch
// 3.00 - Change VVIObjF so it only reads kappa
//
//
//
// Created by Morris Swartz on 7/13/17.
//
// Simplification of VVIObjF by Tamas Vami
//

#ifdef SI_PIXEL_TEMPLATE_STANDALONE
Expand Down Expand Up @@ -696,9 +697,8 @@ int SiPixelTemplateReco2D::PixelTempReco2D(int id,
float sigmaQ = templ2D.sigmavav();
float kappa = templ2D.kappavav();
float xvav = (qtotal / qtfrac0[ipass] - mpv) / sigmaQ;
float beta2 = 1.f;
// VVIObj is a private port of CERNLIB VVIDIS
VVIObjF vvidist(kappa, beta2, 1);
VVIObjF vvidist(kappa);
float prvav = vvidist.fcn(xvav);
probQ = 1.f - prvav;
} else {
Expand Down
6 changes: 3 additions & 3 deletions RecoLocalTracker/SiPixelRecHits/src/SiPixelTemplateSplit.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -24,6 +24,7 @@
// Add shape and charge probabilities for the merged cluster hypothesis (v2.23)
// Incorporate VI-like speed improvements (v2.25)
// Improve speed by eliminating the triple index boost::multiarray objects and add speed switch to optimize the algorithm (v2.30)
// Change VVIObjF so it only reads kappa
//

#ifndef SI_PIXEL_TEMPLATE_STANDALONE
Expand Down Expand Up @@ -131,7 +132,7 @@ int SiPixelTemplateSplit::PixelTempSplit(int id,
float originx, originy, bias, maxpix, minmax;
double chi2x, meanx, chi2y, meany, chi2ymin, chi2xmin, chi21max;
double hchi2, hndof, sigmal1, sigmal2, mpv1, mpv2, arg1, arg2, q05, like, loglike1, loglike2;
double prvav, mpv, sigmaQ, kappa, xvav, beta2;
double prvav, mpv, sigmaQ, kappa, xvav;
float ysum[BYSIZE], xsum[BXSIZE], ysort[BYSIZE], xsort[BXSIZE];
float ysig2[BYSIZE], xsig2[BXSIZE];
float yw2[BYSIZE], xw2[BXSIZE], ysw[BYSIZE], xsw[BXSIZE];
Expand Down Expand Up @@ -602,9 +603,8 @@ int SiPixelTemplateSplit::PixelTempSplit(int id,
assert((sigmaQ > 0.) && (mpv > 0.) && (kappa > 0.01) && (kappa < 10.));
#endif
xvav = ((double)qtotal - mpv) / sigmaQ;
beta2 = 1.;
// VVIObj is a private port of CERNLIB VVIDIS
VVIObj vvidist(kappa, beta2, 1);
VVIObj vvidist(kappa);
prvav = vvidist.fcn(xvav);
prob2Q = 1. - prvav;
if (prob2Q < prob2Qmin) {
Expand Down
85 changes: 85 additions & 0 deletions RecoLocalTracker/SiPixelRecHits/src/VVIObjF.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@
// V1.1 - make dzero call both fcns with a switch
// V1.2 - remove inappriate initializers and add methods to return non-zero/normalized region
// V2.0 - restructuring and speed improvements by V. Innocente
// V3.0 - further simplification and speedup by Tamas Vami
//

#ifndef SI_PIXEL_TEMPLATE_STANDALONE
Expand Down Expand Up @@ -38,6 +39,7 @@ namespace VVIObjFDetails {
//! \param mode - (input) set to 0 to calculate the density function and to 1 to calculate the distribution function
// ***************************************************************************************************************************************

// WARNING: if you change this, dont forget to change VVIObjF::VVIObjF(float kappa) too
VVIObjF::VVIObjF(float kappa, float beta2, int mode) : 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};
Expand Down Expand Up @@ -132,6 +134,89 @@ VVIObjF::VVIObjF(float kappa, float beta2, int mode) : mode_(mode) {

} // VVIObjF

// ***************************************************************************************************************************************
//! Constructor
//! Set Vavilov parameter kappa and calculate the distribution fcn
//! \param kappa - (input) Vavilov kappa parameter [0.01 (Landau-like) < kappa < 10. (Gaussian-like)]
// ***************************************************************************************************************************************

// WARNING: if you change this, dont forget to change the full constructor too
VVIObjF::VVIObjF(float kappa) : mode_(1) {
const float xp[9] = {9.29, 2.47, .89, .36, .15, .07, .03, .02, 0.0};
tvami marked this conversation as resolved.
Show resolved Hide resolved
const float xq[7] = {.012, .03, .08, .26, .87, 3.83, 11.0};
float h_[5];
float 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.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.;
for (k = 1; k < n; ++k) {
l = n - k;
x = omega_ * k;
x1 = h6 * x;
VVIObjFDetails::sincosint(x1, c2, c1);
c1 = vdt::fast_logf(x) - c1;
vdt::fast_sincosf(x1, c3, c4);
xf1 = kappa * (c1 - c4) - x * c2;
xf2 = x * c1 + kappa * (c3 + c2) + t0_ * x;
float s, c;
vdt::fast_sincosf(xf2, s, c);
d1 = q * d * vdt::fast_expf(xf1) / k;
a_[l - 1] = d1 * s;
b_[l - 1] = d1 * c;
a_[n - 1] += q2 * a_[l - 1];
q = -q;
q2 = -q2;
}

} // VVIObjF with kappa only

// *************************************************************************************************************************************
//! Vavilov function method
//! Returns density fcn (mode=0) or distribution fcn (mode=1)
Expand Down