diff --git a/RecoLocalTracker/SiPixelRecHits/interface/VVIObjF.h b/RecoLocalTracker/SiPixelRecHits/interface/VVIObjF.h index 3cc2d89cc0e40..308f67878b22f 100644 --- a/RecoLocalTracker/SiPixelRecHits/interface/VVIObjF.h +++ b/RecoLocalTracker/SiPixelRecHits/interface/VVIObjF.h @@ -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 @@ -23,7 +24,8 @@ // *********************************************************************************************************************** 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, @@ -31,7 +33,6 @@ class VVIObjF { 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_; diff --git a/RecoLocalTracker/SiPixelRecHits/src/SiPixelTemplateReco.cc b/RecoLocalTracker/SiPixelRecHits/src/SiPixelTemplateReco.cc index ac88d5ad51c6f..78d0a97ee20de 100644 --- a/RecoLocalTracker/SiPixelRecHits/src/SiPixelTemplateReco.cc +++ b/RecoLocalTracker/SiPixelRecHits/src/SiPixelTemplateReco.cc @@ -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 @@ -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 diff --git a/RecoLocalTracker/SiPixelRecHits/src/SiPixelTemplateReco2D.cc b/RecoLocalTracker/SiPixelRecHits/src/SiPixelTemplateReco2D.cc index b8cf90678b2ce..35c6c18982ae4 100644 --- a/RecoLocalTracker/SiPixelRecHits/src/SiPixelTemplateReco2D.cc +++ b/RecoLocalTracker/SiPixelRecHits/src/SiPixelTemplateReco2D.cc @@ -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 @@ -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 { diff --git a/RecoLocalTracker/SiPixelRecHits/src/SiPixelTemplateSplit.cc b/RecoLocalTracker/SiPixelRecHits/src/SiPixelTemplateSplit.cc index 80738bae916f2..c67a7aaee0c80 100644 --- a/RecoLocalTracker/SiPixelRecHits/src/SiPixelTemplateSplit.cc +++ b/RecoLocalTracker/SiPixelRecHits/src/SiPixelTemplateSplit.cc @@ -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 @@ -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]; @@ -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) { diff --git a/RecoLocalTracker/SiPixelRecHits/src/VVIObjF.cc b/RecoLocalTracker/SiPixelRecHits/src/VVIObjF.cc index 358a8a22b23b7..a853779235e71 100644 --- a/RecoLocalTracker/SiPixelRecHits/src/VVIObjF.cc +++ b/RecoLocalTracker/SiPixelRecHits/src/VVIObjF.cc @@ -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 @@ -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}; @@ -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}; + 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)