Replace TFormula with FormulaEvaluator in SimpleJectCorrector

Threaded performance measurements using Intel's VTune showed that TFormula use was requiring many locks and was causing significant inefficiencies in CMSSW_7_4 with ROOT 6.02. Although much of the problem goes away in ROOT 6.04 with the new TFormula implementation, moving away from TFormula and to the completely thread efficient reco::FormulaEvaluator also makes construction of the formula thread efficient.
cms-sw · Oct 1, 2015 · 83aa621 · 83aa621
1 parent 70fc7bd
commit 83aa621
Show file tree

Hide file tree

Showing 3 changed files with 14 additions and 80 deletions.
diff --git a/CondFormats/JetMETObjects/BuildFile.xml b/CondFormats/JetMETObjects/BuildFile.xml
@@ -1,6 +1,7 @@
 <use   name="DataFormats/CaloTowers"/>
 <use   name="FWCore/ParameterSet"/>
 <use   name="FWCore/Utilities"/>
+<use   name="CommonTools/Utils"/>
 <use   name="CondFormats/Serialization"/>
 <use   name="CondFormats/DataRecord"/>
 <use   name="boost_serialization"/>

diff --git a/CondFormats/JetMETObjects/interface/SimpleJetCorrector.h b/CondFormats/JetMETObjects/interface/SimpleJetCorrector.h
@@ -6,21 +6,18 @@
 #include <string>
 #include <vector>
 
-#include <TFormula.h>
-#include <RVersion.h>
 #include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
 
+#include "CommonTools/Utils/interface/FormulaEvaluator.h"
+
 class JetCorrectorParameters;
 
 class SimpleJetCorrector 
 {
  public:
   //-------- Constructors --------------
-  SimpleJetCorrector();
   SimpleJetCorrector(const std::string& fDataFile, const std::string& fOption = "");
   SimpleJetCorrector(const JetCorrectorParameters& fParameters);
-  //-------- Destructor -----------------
-  ~SimpleJetCorrector();
   //-------- Member functions -----------
   void   setInterpolation(bool fInterpolation) {mDoInterpolation = fInterpolation;}
   float  correction(const std::vector<float>& fX,const std::vector<float>& fY) const;  
@@ -30,17 +27,13 @@ class SimpleJetCorrector
   //-------- Member functions -----------
   SimpleJetCorrector(const SimpleJetCorrector&);
   SimpleJetCorrector& operator= (const SimpleJetCorrector&);
-#if ROOT_VERSION_CODE >= ROOT_VERSION(6,03,00)
-  float    invert(const Double_t *args, const Double_t *params) const;
-#else
-  float    invert(const std::vector<float>& fX, TFormula&) const;
-#endif
+  float    invert(const double *args, const double *params) const;
   float    correctionBin(unsigned fBin,const std::vector<float>& fY) const;
   unsigned findInvertVar();
   void     setFuncParameters();
   //-------- Member variables -----------
   JetCorrectorParameters  mParameters;
-  TFormula                mFunc;
+  reco::FormulaEvaluator  mFunc;
   unsigned                mInvertVar; 
   bool                    mDoInterpolation;
 };

diff --git a/CondFormats/JetMETObjects/src/SimpleJetCorrector.cc b/CondFormats/JetMETObjects/src/SimpleJetCorrector.cc
@@ -5,21 +5,13 @@
 #include <sstream>
 #include <cmath>
 
-//------------------------------------------------------------------------
-//--- Default SimpleJetCorrector constructor -----------------------------
-//------------------------------------------------------------------------
-SimpleJetCorrector::SimpleJetCorrector()
-{
-  mDoInterpolation = false;
-  mInvertVar       = 9999;
-}
 //------------------------------------------------------------------------
 //--- SimpleJetCorrector constructor -------------------------------------
 //--- reads arguments from a file ----------------------------------------
 //------------------------------------------------------------------------
 SimpleJetCorrector::SimpleJetCorrector(const std::string& fDataFile, const std::string& fOption):
   mParameters(fDataFile,fOption),
-  mFunc("function",((mParameters.definitions()).formula()).c_str())
+  mFunc((mParameters.definitions()).formula())
 {
   mDoInterpolation = false;
   if (mParameters.definitions().isResponse())
@@ -31,18 +23,12 @@ SimpleJetCorrector::SimpleJetCorrector(const std::string& fDataFile, const std::
 //------------------------------------------------------------------------
 SimpleJetCorrector::SimpleJetCorrector(const JetCorrectorParameters& fParameters):
   mParameters(fParameters),
-  mFunc("function",((mParameters.definitions()).formula()).c_str())
+  mFunc((mParameters.definitions()).formula())
 {
   mDoInterpolation = false;
   if (mParameters.definitions().isResponse())
     mInvertVar = findInvertVar();
 }
-//------------------------------------------------------------------------
-//--- SimpleJetCorrector destructor --------------------------------------
-//------------------------------------------------------------------------
-SimpleJetCorrector::~SimpleJetCorrector()
-{
-}
 
 //------------------------------------------------------------------------
 //--- calculates the correction ------------------------------------------
@@ -106,42 +92,20 @@ float SimpleJetCorrector::correctionBin(unsigned fBin,const std::vector<float>&
     }
   const std::vector<float>& par = mParameters.record(fBin).parameters();
 
-#if ROOT_VERSION_CODE >= ROOT_VERSION(6,03,00)
-  Double_t params[par.size() - 2 * N];
+  double params[par.size() - 2 * N];
   for(unsigned int i=2*N;i<par.size();i++)
     {
       params[i-2*N] = par[i];
     }
-  Double_t x[4] = {};
+  double x[4] = {};
   for(unsigned i=0;i<N;i++)
     {
       x[i] = (fY[i] < par[2*i]) ? par[2*i] : (fY[i] > par[2*i+1]) ? par[2*i+1] : fY[i];
     }
   if (mParameters.definitions().isResponse()) {
     return invert(x, params);
   }
-  return mFunc.EvalPar(x, params);
-#else
-  float result = -1;
-  //Have to do calculation using a temporary TFormula to avoid
-  // thread safety issues
-  TFormula tFunc(mFunc);
-
-  for(unsigned int i=2*N;i<par.size();i++)
-    tFunc.SetParameter(i-2*N,par[i]);
-  float x[4] = {};
-  std::vector<float> tmp;
-  for(unsigned i=0;i<N;i++)
-    {
-      x[i] = (fY[i] < par[2*i]) ? par[2*i] : (fY[i] > par[2*i+1]) ? par[2*i+1] : fY[i];
-      tmp.push_back(x[i]);
-    }
-  if (mParameters.definitions().isResponse())
-    result = invert(tmp,tFunc);
-  else
-    result = tFunc.Eval(x[0],x[1],x[2],x[3]);
-  return result;
-#endif
+  return mFunc.evaluate(reco::formula::ArrayAdaptor(x,N), reco::formula::ArrayAdaptor(params,par.size()-2*N) );
 }
 //------------------------------------------------------------------------
 //--- find invertion variable (JetPt) ------------------------------------
@@ -163,49 +127,25 @@ unsigned SimpleJetCorrector::findInvertVar()
 //------------------------------------------------------------------------
 //--- inversion ----------------------------------------------------------
 //------------------------------------------------------------------------
-#if ROOT_VERSION_CODE >= ROOT_VERSION(6,03,00)
-float SimpleJetCorrector::invert(const Double_t *args, const Double_t *params) const
+float SimpleJetCorrector::invert(const double *args, const double *params) const
 {
   unsigned nMax = 50;
   float precision = 0.0001;
   float rsp = 1.0;
   float e = 1.0;
-  Double_t x[4];
+  double x[4];
   unsigned nLoop=0;
 
   // 4 dimensions (x, y, z, t) supported in TFormula
-  memcpy(&x, args, sizeof(Double_t) * 4);
+  memcpy(&x, args, sizeof(double) * 4);
 
   while(e > precision && nLoop < nMax)
     {
-      rsp = mFunc.EvalPar(x, params);
+      rsp = mFunc.evaluate(reco::formula::ArrayAdaptor(x,4), reco::formula::ArrayAdaptor(params,mFunc.numberOfParameters()));
       float tmp = x[mInvertVar] * rsp;
       e = fabs(tmp - args[mInvertVar])/args[mInvertVar];
       x[mInvertVar] = args[mInvertVar]/rsp;
       nLoop++;
     }
   return 1./rsp;
 }
-#else
-float SimpleJetCorrector::invert(const std::vector<float>& fX, TFormula& tFunc) const
-{
-  unsigned nMax = 50;
-  unsigned N = fX.size();
-  float precision = 0.0001;
-  float rsp = 1.0;
-  float e = 1.0;
-  float x[4] = {0.0,0.0,0.0,0.0};
-  for(unsigned i=0;i<N;i++)
-    x[i] = fX[i];
-  unsigned nLoop=0;
-  while(e > precision && nLoop < nMax)
-    {
-      rsp = tFunc.Eval(x[0],x[1],x[2],x[3]);
-      float tmp = x[mInvertVar] * rsp;
-      e = fabs(tmp - fX[mInvertVar])/fX[mInvertVar];
-      x[mInvertVar] = fX[mInvertVar]/rsp;
-      nLoop++;
-    }
-  return 1./rsp;
-}
-#endif