CMS 3D CMS Logo

SimpleJetCorrector.cc
Go to the documentation of this file.
4 #include <iostream>
5 #include <sstream>
6 #include <cmath>
7 
8 //------------------------------------------------------------------------
9 //--- SimpleJetCorrector constructor -------------------------------------
10 //--- reads arguments from a file ----------------------------------------
11 //------------------------------------------------------------------------
13  mParameters(fDataFile,fOption),
14  mFunc((mParameters.definitions()).formula())
15 {
16  mDoInterpolation = false;
19 }
20 //------------------------------------------------------------------------
21 //--- SimpleJetCorrector constructor -------------------------------------
22 //--- reads arguments from a file ----------------------------------------
23 //------------------------------------------------------------------------
25  mParameters(fParameters),
27 {
28  mDoInterpolation = false;
31 }
32 
33 //------------------------------------------------------------------------
34 //--- calculates the correction ------------------------------------------
35 //------------------------------------------------------------------------
36 float SimpleJetCorrector::correction(const std::vector<float>& fX,const std::vector<float>& fY) const
37 {
38  float result = 1.;
39  float tmp = 0.0;
40  float cor = 0.0;
41  int bin = -1;
42  bin = (fX.size()<=3 && !fX.empty()) ? mParameters.binIndexN(fX) : mParameters.binIndex(fX);
43  if (bin<0)
44  return result;
45  if (!mDoInterpolation)
46  result = correctionBin(bin,fY);
47  else
48  {
49  for(unsigned i=0;i<mParameters.definitions().nBinVar();i++)
50  {
51  float xMiddle[3];
52  float xValue[3];
53  int prevBin = mParameters.neighbourBin((unsigned)bin,i,false);
54  int nextBin = mParameters.neighbourBin((unsigned)bin,i,true);
55  if (prevBin>=0 && nextBin>=0)
56  {
57  xMiddle[0] = mParameters.record(prevBin).xMiddle(i);
58  xMiddle[1] = mParameters.record(bin).xMiddle(i);
59  xMiddle[2] = mParameters.record(nextBin).xMiddle(i);
60  xValue[0] = correctionBin(prevBin,fY);
61  xValue[1] = correctionBin(bin,fY);
62  xValue[2] = correctionBin(nextBin,fY);
63  cor = quadraticInterpolation(fX[i],xMiddle,xValue);
64  tmp+=cor;
65  }
66  else
67  {
68  cor = correctionBin(bin,fY);
69  tmp+=cor;
70  }
71  }
72  result = tmp/mParameters.definitions().nBinVar();
73  }
74  return result;
75 }
76 //------------------------------------------------------------------------
77 //--- calculates the correction for a specific bin -----------------------
78 //------------------------------------------------------------------------
79 float SimpleJetCorrector::correctionBin(unsigned fBin,const std::vector<float>& fY) const
80 {
81  if (fBin >= mParameters.size())
82  {
83  std::stringstream sserr;
84  sserr<<"wrong bin: "<<fBin<<": only "<<mParameters.size()<<" available!";
85  handleError("SimpleJetCorrector",sserr.str());
86  }
87  unsigned N = fY.size();
88  if (N > 4)
89  {
90  std::stringstream sserr;
91  sserr<<"two many variables: "<<N<<" maximum is 4";
92  handleError("SimpleJetCorrector",sserr.str());
93  }
94  const std::vector<float>& par = mParameters.record(fBin).parameters();
95 
96  double params[par.size() - 2 * N];
97  for(unsigned int i=2*N;i<par.size();i++)
98  {
99  params[i-2*N] = par[i];
100  }
101  double x[4] = {};
102  for(unsigned i=0;i<N;i++)
103  {
104  x[i] = (fY[i] < par[2*i]) ? par[2*i] : (fY[i] > par[2*i+1]) ? par[2*i+1] : fY[i];
105  }
107  return invert(x, params);
108  }
109  return mFunc.evaluate(reco::formula::ArrayAdaptor(x,N), reco::formula::ArrayAdaptor(params,par.size()-2*N) );
110 }
111 //------------------------------------------------------------------------
112 //--- find invertion variable (JetPt) ------------------------------------
113 //------------------------------------------------------------------------
115 {
116  unsigned result = 9999;
117  std::vector<std::string> vv = mParameters.definitions().parVar();
118  for(unsigned i=0;i<vv.size();i++)
119  if (vv[i]=="JetPt")
120  {
121  result = i;
122  break;
123  }
124  if (result >= vv.size())
125  handleError("SimpleJetCorrector","Response inversion is required but JetPt is not specified as parameter");
126  return result;
127 }
128 //------------------------------------------------------------------------
129 //--- inversion ----------------------------------------------------------
130 //------------------------------------------------------------------------
131 float SimpleJetCorrector::invert(const double *args, const double *params) const
132 {
133  unsigned nMax = 50;
134  float precision = 0.0001;
135  float rsp = 1.0;
136  float e = 1.0;
137  double x[4];
138  unsigned nLoop=0;
139 
140  // 4 dimensions (x, y, z, t) supported in TFormula
141  memcpy(&x, args, sizeof(double) * 4);
142 
143  while(e > precision && nLoop < nMax)
144  {
146  float tmp = x[mInvertVar] * rsp;
147  e = fabs(tmp - args[mInvertVar])/args[mInvertVar];
148  x[mInvertVar] = args[mInvertVar]/rsp;
149  nLoop++;
150  }
151  return 1./rsp;
152 }
unsigned int numberOfParameters() const
float correctionBin(unsigned fBin, const std::vector< float > &fY) const
float invert(const double *args, const double *params) const
std::vector< float > parameters() const
const Definitions & definitions() const
SimpleJetCorrector(const std::string &fDataFile, const std::string &fOption="")
const Record & record(unsigned fBin) const
std::vector< std::string > parVar() const
int binIndex(const std::vector< float > &fX) const
double evaluate(V const &iVariables, P const &iParameters) const
float xMiddle(unsigned fVar) const
bin
set the eta bin as selection string.
#define N
Definition: blowfish.cc:9
int binIndexN(const std::vector< float > &fX) const
std::vector< std::vector< double > > tmp
Definition: MVATrainer.cc:100
float correction(const std::vector< float > &fX, const std::vector< float > &fY) const
reco::FormulaEvaluator mFunc
int neighbourBin(unsigned fIndex, unsigned fVar, bool fNext) const
JetCorrectorParameters mParameters