HcalAlgos/interface/genlkupmap.h

0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002 #include "CalibCalorimetry/HcalAlgos/interface/HcalShapeIntegrator.h"
0003 #include <algorithm>  // for "max","min"
0004 #include <cmath>
0005 #include <iostream>
0006 #include <memory>
0007
0008 // Function generates a lookup map for a passed-in function (via templated object algoObject,
0009 // which must contain method "calcpair" that spits out (x,y) pair from a type float seed.
0010 // Each map y-value is separated from the previous value by a programmable fractional error
0011 // relative to the previous value.
0012 //
0013 // Currently this function coded for only monotonically increasing or
0014 // decreasing functions...
0015
0016 template <class T_algo>
0017 void genlkupmap(double smin,
0018                 double smax,
0019                 double max_fracerror,
0020                 double min_sstep,
0021                 T_algo& algoObject,
0022                 std::map<double, double>& m_ylookup) {
0023   std::pair<double, double> thisxy, lastxy, laststoredxy;
0024   std::pair<double, double> minxy = algoObject.calcpair(smin);
0025   std::pair<double, double> maxxy = algoObject.calcpair(smax);
0026
0027   double slope = maxxy.second - minxy.second;
0028   slope = (slope >= 0.0) ? 1.0 : -1.0;
0029
0030   double sstep = min_sstep;
0031
0032   for (double s = smin; lastxy.first < smax; s += sstep) {
0033     thisxy = algoObject.calcpair(s);
0034
0035     double fracerror = slope * (thisxy.second - laststoredxy.second) / thisxy.second;
0036     double fracchange = slope * (thisxy.second - lastxy.second) / thisxy.second;
0037
0038     bool store_cur_pair = false;
0039     bool store_prev_pair = false;
0040
0041 #if 0
0042     char str[80];
0043     sprintf(str, "%7.1f %7.1f (%8.3f %8.4f) %8.5f %8.5f",
0044         s, sstep, thisxy.first, thisxy.second, fracerror, fracchange);
0045     cout << str;
0046 #endif
0047
0048     if (s == smin) {
0049       store_cur_pair = true;
0050     } else if ((fracerror > 2 * max_fracerror) || (thisxy.first > smax)) {
0051       if (sstep > min_sstep) {
0052         // possibly overstepped the next entry, back up and reduce the step size
0053         s -= sstep;
0054         sstep = std::max(0.5 * sstep, min_sstep);
0055 #if 0
0056     cout << endl;
0057 #endif
0058         continue;
0059       } else if (lastxy.second == laststoredxy.second) {
0060         store_cur_pair = true;
0061
0062         // current step size is too big for the max allowed fractional error,
0063         // store current value and issue a warning.
0064         //
0065         //  edm::LogWarning("HcalPulseContainmentCorrection::genlkupmap") << " fractional error max exceeded";
0066       } else {
0067         store_prev_pair = true;
0068
0069         // re-evaluate current yval with prev yval.
0070         fracerror = slope * (thisxy.second - lastxy.second) / thisxy.second;
0071
0072         if (fracerror > 2 * max_fracerror) {
0073           store_cur_pair = true;
0074
0075           // current step size is too big for the max allowed fractional error,
0076           // store current value and issue a warning.
0077           //
0078           //      edm::LogWarning("HcalPulseContainmentCorrection::genlkupmap") << " fractional error max exceeded";
0079         }
0080       }
0081     } else if ((fracerror < 1.9 * max_fracerror) && (fracchange < 0.1 * max_fracerror) &&
0082                (thisxy.first < 0.99 * smax)) {
0083       // adapt step size to reduce iterations
0084       sstep *= 2.0;
0085     }
0086
0087     if (thisxy.first > smax)
0088       store_cur_pair = true;
0089
0090     if (store_prev_pair) {
0091       m_ylookup[lastxy.first] = lastxy.second;
0092       laststoredxy = lastxy;
0093     }
0094     if (store_cur_pair) {
0095       m_ylookup[thisxy.first] = thisxy.second;
0096       laststoredxy = thisxy;
0097     }
0098
0099     lastxy = thisxy;
0100
0101 #if 0
0102     sprintf(str, " %c %c",
0103         store_cur_pair ? 'C' : ' ',
0104         store_prev_pair ? 'P' : ' ');
0105     cout << str << endl;
0106 #endif
0107   }
0108 }
0109
0110 //======================================================================
0111 // Fixed Phase mode amplitude correction factor generation routines:
0112
0113 template <class S>
0114 class RecoFCcorFactorAlgo {
0115 public:
0116   RecoFCcorFactorAlgo(int num_samples, double fixedphase_ns);
0117   std::pair<double, double> calcpair(double);
0118
0119 private:
0120   double fixedphasens_;
0121   double integrationwindowns_;
0122   double time0shiftns_;
0123   S shape_;
0124   const std::unique_ptr<HcalShapeIntegrator> integrator_;
0125 };