epic-lfhcal-tbana/NewStructure/AdaptiveLangau.h

0001 #ifndef LFHCAL_ADAPTIVE_LANGAU_H
0002 #define LFHCAL_ADAPTIVE_LANGAU_H
0003
0004 // Numerical helper only: no minimizer, fit limits, selection, or ROOT objects.
0005 // Extracted from LangauNumerics.h at 3a643646eb301c7112449e8f6d6ab31f12c72943.
0006 // The subdivision rule, breakpoints, tolerances, and +/-5 sigma convention
0007 // are retained. Expensive cross-checks belong in tests, not every production fit.
0008 #include <algorithm>
0009 #include <array>
0010 #include <cmath>
0011 #include <functional>
0012 #include <stdexcept>
0013 #include <vector>
0014
0015 namespace lfhcal { namespace adaptive_langau {
0016 using Function = std::function<double(double)>;
0017 struct Parameters { double width, mp, area, sigma; };
0018 struct IntegrationConfig {
0019   double relative_tolerance = 1e-10;
0020   double gaussian_span = 5.;
0021   int max_calls = 500000;
0022   int max_depth = 18;
0023   double absolute_tolerance_fraction = .01;
0024 };
0025 struct PeakWidth { double peak, left, right, fwhm, height; };
0026
0027 inline void require(bool ok, const char* message) {
0028   if (!ok) throw std::runtime_error(message);
0029 }
0030 inline double checked(const Function& f, double x) {
0031   const double y = f(x);
0032   require(std::isfinite(y) && y >= 0., "nonfinite or negative Langau evaluation");
0033   return y;
0034 }
0035
0036 namespace detail {
0037 // Positive nodes/weights of the 16-point Gauss-Legendre rule on [-1,1].
0038 constexpr std::array<double,8> nodes = {
0039   .095012509837637440185, .28160355077925891323,
0040   .45801677765722738634, .61787624440264374845,
0041   .75540440835500303390, .86563120238783174388,
0042   .94457502307323257608, .98940093499164993260};
0043 constexpr std::array<double,8> weights = {
0044   .18945061045506849629, .18260341504492358887,
0045   .16915651939500253819, .14959598881657673208,
0046   .12462897125553387205, .095158511682492784810,
0047   .062253523938647892863, .027152459411754094852};
0048 constexpr std::array<double,16> breakpoints = {
0049   -16.,-8.,-5.5,-4.,-2.,-1.,0.,1.,2.,4.,5.,8.,12.,16.,50.,300.};
0050
0051 inline double gauss16(const Function& f, double a, double b,
0052                       int& calls, int maximum) {
0053   calls += 16;
0054   require(calls <= maximum, "Langau quadrature evaluation budget exhausted");
0055   const double c=(a+b)/2, d=(b-a)/2;
0056   double sum=0.;
0057   for (std::size_t i=0; i<nodes.size(); ++i)
0058     sum += weights[i]*(checked(f,c-d*nodes[i])+checked(f,c+d*nodes[i]));
0059   return d*sum;
0060 }
0061 inline double subdivide(const Function& f, double a, double b, double coarse,
0062                          double atol, double rtol, int depth,
0063                          int& calls, int maximum) {
0064   const double mid=(a+b)/2;
0065   const double l=gauss16(f,a,mid,calls,maximum);
0066   const double r=gauss16(f,mid,b,calls,maximum);
0067   const double fine=l+r;
0068   if (std::abs(fine-coarse) <= atol+rtol*std::abs(fine)) return fine;
0069   require(depth>0 && a<mid && mid<b, "Langau quadrature did not converge");
0070   return subdivide(f,a,mid,l,atol/2,rtol,depth-1,calls,maximum)
0071        + subdivide(f,mid,b,r,atol/2,rtol,depth-1,calls,maximum);
0072 }
0073 } // namespace detail
0074
0075 // u=(t-mpc)/width resolves the narrow Landau on its own scale. A single
0076 // adaptive integral over the unsplit interval could miss its entire core.
0077 class Convolution {
0078  public:
0079   Convolution(Parameters p, Function standardLandau, IntegrationConfig cfg = {})
0080       : p_(p), pdf_(standardLandau), cfg_(cfg) {
0081     require(std::isfinite(p.width) && p.width>0 && std::isfinite(p.sigma) && p.sigma>0
0082             && std::isfinite(p.mp) && std::isfinite(p.area) && p.area>0,
0083             "invalid Langau parameters");
0084     require(cfg.relative_tolerance>0 && cfg.relative_tolerance<1
0085             && cfg.gaussian_span>=5 && cfg.gaussian_span<=12
0086             && cfg.max_calls>0 && cfg.max_depth>0
0087             && cfg.absolute_tolerance_fraction>0, "invalid Langau integration settings");
0088   }
0089   double operator()(double x) const {
0090     require(std::isfinite(x), "nonfinite Langau coordinate");
0091     const double mpc=p_.mp+.22278298*p_.width;
0092     const double lo=(x-cfg_.gaussian_span*p_.sigma-mpc)/p_.width;
0093     const double hi=(x+cfg_.gaussian_span*p_.sigma-mpc)/p_.width;
0094     std::vector<double> edges{lo,hi};
0095     const auto add=[&](double u) { if (u>lo && u<hi) edges.push_back(u); };
0096     for (double u:detail::breakpoints) add(u);
0097     const double extent=std::max(std::abs(lo),std::abs(hi));
0098     require(std::isfinite(extent) && extent<1e12, "Langau parameter ratio exceeds numerical domain");
0099     for (double u=32.; u<extent; u*=2.) { add(u); add(-u); }
0100     for (int k=-static_cast<int>(cfg_.gaussian_span);
0101          k<=static_cast<int>(cfg_.gaussian_span); ++k)
0102       add((x+k*p_.sigma-mpc)/p_.width);
0103     std::sort(edges.begin(),edges.end());
0104     edges.erase(std::unique(edges.begin(),edges.end()),edges.end());
0105     Function integrand=[&](double u) {
0106       const double z=(x-mpc-p_.width*u)/p_.sigma;
0107       return pdf_(u)*std::exp(-.5*z*z);
0108     };
0109     int calls=0;
0110     double sum=0.;
0111     const double atol=cfg_.relative_tolerance*cfg_.absolute_tolerance_fraction/(edges.size()-1);
0112     for (std::size_t i=1; i<edges.size(); ++i) {
0113       const double a=edges[i-1], b=edges[i];
0114       const double coarse=detail::gauss16(integrand,a,b,calls,cfg_.max_calls);
0115       sum += detail::subdivide(integrand,a,b,coarse,atol,cfg_.relative_tolerance,
0116                                cfg_.max_depth,calls,cfg_.max_calls);
0117     }
0118     const double value=p_.area*sum/(p_.sigma*std::sqrt(2*std::acos(-1.)));
0119     require(std::isfinite(value) && value>=0., "invalid Langau integral");
0120     return value;
0121   }
0122  private:
0123   Parameters p_;
0124   Function pdf_;
0125   IntegrationConfig cfg_;
0126 };
0127
0128 // Same bounded peak/half-height procedure as the tested adaptive helper.
0129 // It evaluates the continuous model, not a sampled/FFT graph. This is numerical
0130 // extraction of ScaleH/FWHM, not an additional fit or a goodness-of-fit cut.
0131 inline PeakWidth peakWidth(const Function& f, double centre, double scale) {
0132   require(std::isfinite(centre) && std::isfinite(scale) && scale>0.,
0133           "invalid Langau peak-search scale");
0134   constexpr int scanBins=160;
0135   double a=centre-8*scale, b=centre+12*scale, lo=0., hi=0.;
0136   bool bracketed=false;
0137   for (int expand=0; expand<8; ++expand) {
0138     std::vector<double> y(scanBins+1);
0139     int best=0;
0140     for (int i=0; i<=scanBins; ++i) {
0141       y[i]=checked(f,a+(b-a)*i/scanBins);
0142       if (y[i]>y[best]) best=i;
0143     }
0144     require(y[best]>0., "zero curve across Langau peak-search window");
0145     if (best>0 && best<scanBins) {
0146       int peaks=0;
0147       for (int i=1; i<scanBins; ++i)
0148         if (y[i]>y[i-1] && y[i]>=y[i+1] && y[i]>1e-5*y[best]) ++peaks;
0149       require(peaks==1, "multiple resolved peaks in Langau model");
0150       lo=a+(b-a)*(best-1)/scanBins;
0151       hi=a+(b-a)*(best+1)/scanBins;
0152       bracketed=true;
0153       break;
0154     }
0155     a=centre+2*(a-centre); b=centre+2*(b-centre);
0156   }
0157   require(bracketed, "unable to bracket Langau peak");
0158   const double golden=(std::sqrt(5.)-1)/2, xtol=1e-8*std::max(1.,scale);
0159   double l=hi-golden*(hi-lo), r=lo+golden*(hi-lo);
0160   double fl=checked(f,l), fr=checked(f,r);
0161   for (int k=0; k<160 && hi-lo>xtol; ++k) {
0162     if (fl>fr) {hi=r; r=l; fr=fl; l=hi-golden*(hi-lo); fl=checked(f,l);}
0163     else {lo=l; l=r; fl=fr; r=lo+golden*(hi-lo); fr=checked(f,r);}
0164   }
0165   require(hi-lo<=xtol, "Langau peak search did not converge");
0166   const double peak=(lo+hi)/2, height=checked(f,peak);
0167   require(height>0., "nonpositive Langau peak height");
0168   const auto crossing=[&](int direction) {
0169     double inside=peak, outside=peak, distance=scale;
0170     bool found=false;
0171     for (int k=0; k<64; ++k) {
0172       outside=peak+direction*distance;
0173       if (checked(f,outside)<height/2) {found=true; break;}
0174       distance*=1.8;
0175     }
0176     require(found, "unable to bracket Langau half height");
0177     for (int k=0; k<200; ++k) {
0178       const double mid=(inside+outside)/2, residual=checked(f,mid)/height-.5;
0179       if (std::abs(residual)<1e-10) return mid;
0180       if (residual>0) inside=mid; else outside=mid;
0181       if (std::abs(inside-outside)<xtol*.01) return (inside+outside)/2;
0182     }
0183     throw std::runtime_error("Langau half-height search did not converge");
0184   };
0185   const double left=crossing(-1), right=crossing(1);
0186   require(left<peak && peak<right && std::abs(checked(f,left)/height-.5)<1e-8
0187           && std::abs(checked(f,right)/height-.5)<1e-8, "invalid Langau half-height crossings");
0188   return {peak,left,right,right-left,height};
0189 }
0190 }} // namespace lfhcal::adaptive_langau
0191 #endif