genlkupmap.h File Reference

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include <algorithm>
#include <math.h>
#include <iostream>

Go to the source code of this file.

Classes
class	RecoFCcorFactorAlgo< S >

Functions
template<class T_algo >
void	genlkupmap (double smin, double smax, double max_fracerror, double min_sstep, T_algo &algoObject, std::map< double, double > &m_ylookup)

Function Documentation

template<class T_algo >

void genlkupmap	(	double	smin,
		double	smax,
		double	max_fracerror,
		double	min_sstep,
		T_algo &	algoObject,
		std::map< double, double > &	m_ylookup
	)

Definition at line 15 of file genlkupmap.h.

References gather_cfg::cout, if(), max(), asciidump::s, and slope.

{
  std::pair<double,double> thisxy, lastxy, laststoredxy;
  std::pair<double,double> minxy = algoObject.calcpair(smin);
  std::pair<double,double> maxxy = algoObject.calcpair(smax);
  
  double slope = maxxy.second - minxy.second;
  slope = (slope >= 0.0) ? 1.0 : -1.0;

  double sstep = min_sstep;

  for (double s=smin; lastxy.first<smax; s += sstep) {

    thisxy = algoObject.calcpair(s);

    double fracerror  = slope*(thisxy.second - laststoredxy.second)/ thisxy.second;
    double fracchange = slope*(thisxy.second - lastxy.second)/ thisxy.second;

    bool store_cur_pair  = false;
    bool store_prev_pair = false;

#if 0
    char str[80];
    sprintf(str, "%7.1f %7.1f (%8.3f %8.4f) %8.5f %8.5f",
        s, sstep, thisxy.first, thisxy.second, fracerror, fracchange);
    cout << str;
#endif

    if (s == smin) {
      store_cur_pair = true;
    }
    else if ((fracerror > 2*max_fracerror) ||
             (thisxy.first > smax) ) {
      if (sstep > min_sstep) {
        // possibly overstepped the next entry, back up and reduce the step size
        s -= sstep;
        sstep = std::max(0.5*sstep, min_sstep);
#if 0
    cout << endl;
#endif
        continue;
      }
      else if (lastxy.second == laststoredxy.second) {
        store_cur_pair = true;

        // current step size is too big for the max allowed fractional error,
        // store current value and issue a warning.
        //
    edm::LogWarning("HcalPulseContainmentCorrection::genlkupmap") << " fractional error max exceeded";
      }
      else {
        store_prev_pair = true;

        // re-evaluate current yval with prev yval.
        fracerror = slope*(thisxy.second - lastxy.second)/ thisxy.second;

        if (fracerror > 2*max_fracerror) {
          store_cur_pair = true;

          // current step size is too big for the max allowed fractional error,
          // store current value and issue a warning.
          //
      edm::LogWarning("HcalPulseContainmentCorrection::genlkupmap") << " fractional error max exceeded";
        }
      }
    }
    else if ((fracerror  < 1.9*max_fracerror) &&
             (fracchange < 0.1*max_fracerror) &&
             (thisxy.first < 0.99 * smax) ) {
      // adapt step size to reduce iterations
      sstep *= 2.0;
    }

    if (thisxy.first > smax)
      store_cur_pair = true;
      
    if (store_prev_pair) {
      m_ylookup[lastxy.first] = lastxy.second;
      laststoredxy            = lastxy;
    }
    if (store_cur_pair) {
      m_ylookup[thisxy.first] = thisxy.second;
      laststoredxy            = thisxy;
    }

    lastxy = thisxy;
    
#if 0
    sprintf(str, " %c %c",
        store_cur_pair ? 'C' : ' ',
        store_prev_pair ? 'P' : ' ');
    cout << str << endl;
#endif

  }

}