CMS 3D CMS Logo

JetCorrectorParametersHelper.cc
Go to the documentation of this file.
1 //
2 // Original Author: Alexx Perloff Feb 22, 2017
3 // $Id: JetCorrectorParameters.cc,v 1.20 2012/03/01 18:24:53 srappocc Exp $
4 //
5 // Helper class for JetCorrectorParameters
6 //
8 #include <sstream>
9 #include <cstdlib>
10 #include <algorithm>
11 #include <cmath>
12 #include <iterator>
13 
14 //------------------------------------------------------------------------
15 //--- JetCorrectorParameters::JetCorrectorParametersHelper initializer ---
16 //--- initializes the mBinMap for quick lookup of mRecords ---------------
17 //------------------------------------------------------------------------
19 {
20  mIndexMap.clear();
21  mMap.clear();
22  mBinBoundaries.clear();
24 }
26  const std::vector<JetCorrectorParameters::Record>& mRecordsLocal)
27 {
28  SIZE = mDefinitionsLocal.nBinVar();
30  {
31  std::stringstream sserr;
32  sserr<<"The number of binned variables requested ("<<SIZE<<") is greater than the number allowed ("
34  handleError("JetCorrectorParametersHelper",sserr.str());
35  }
36 
38  size_t start=0, end=0;
39  size_t nRec = mRecordsLocal.size();
40  size_t indexMapSize=0, tmpIndexMapSize=0;
41  for (unsigned i = 0; i < nRec; ++i)
42  {
43  for (unsigned j=0;j<SIZE;j++)
44  {
45  if(j<SIZE-1 && std::find(mBinBoundaries[j].begin(),mBinBoundaries[j].end(),mRecordsLocal[i].xMin(j))==mBinBoundaries[j].end())
46  mBinBoundaries[j].push_back(mRecordsLocal[i].xMin(j));
47  else if(j==SIZE-1) {
48  if(i==0)
49  mBinBoundaries[j].reserve(mRecordsLocal.size());
50 
51  mBinBoundaries[j].push_back(mRecordsLocal[i].xMin(j));
52 
53  if(SIZE>1 && (i==nRec-1 || mRecordsLocal[i].xMin(j-1)!=mRecordsLocal[i+1].xMin(j-1))) {
54  end = i;
55  //mMap.emplace(gen_tuple<SIZE-1>([&](size_t k){return mRecordsLocal[i].xMin(k);}),std::make_pair(start,end));
56  mMap.emplace(gen_tuple<JetCorrectorParameters::MAX_SIZE_DIMENSIONALITY-1>([&](size_t k){return (k<SIZE-1)?mRecordsLocal[i].xMin(k):-9999;}),std::make_pair(start,end));
57  start = i+1;
58  }
59  }
60  }
61  indexMapSize = mIndexMap.size();
62  tuple_type tmpTuple = gen_tuple<JetCorrectorParameters::MAX_SIZE_DIMENSIONALITY>([&](size_t k){return (k<SIZE)?mRecordsLocal[i].xMin(k):-9999;});
63  mIndexMap.emplace(tmpTuple,i);
64  tmpIndexMapSize = mIndexMap.size();
65  if(indexMapSize==tmpIndexMapSize)
66  {
67  size_t existing_index = mIndexMap.find(tmpTuple)->second;
68  std::stringstream sserr;
69  sserr<<"Duplicate binning in record found (existing index,current index)=("
70  <<existing_index<<","<<i<<")"<<std::endl<<"\tBins(lower bounds)="<<tmpTuple;
71  handleError("JetCorrectorParametersHelper",sserr.str());
72  }
73  }
74  if (mBinBoundaries[SIZE-1].size()!=nRec)
75  {
76  std::stringstream sserr;
77  sserr<<"Did not find all bin boundaries for dimension "<<SIZE-1<<"!!!"<<std::endl
78  <<"Found "<<mBinBoundaries[SIZE-1].size()<<" out of "<<nRec<<" records";
79  handleError("JetCorrectorParametersHelper",sserr.str());
80  }
81  indexMapSize = mIndexMap.size();
82  if (indexMapSize!=nRec)
83  {
84  handleError("JetCorrectorParametersHelper","The mapping of bin lower bounds to indices does not contain all possible entries!!!");
85  }
86  // This function checks that all of the middle binned parameters (not the first or last) contain the same range of bins.
87  // If not an exception will be thrown as the mapping will work only if all bins but the last are uniform.
88  if (SIZE>2) checkMiddleBinUniformity(mRecordsLocal);
89 }
90 
91 void JetCorrectorParametersHelper::checkMiddleBinUniformity(const std::vector<JetCorrectorParameters::Record>& mRecords) const
92 {
93  unsigned N = SIZE-2;
94  size_t nRec = mRecords.size();
95  std::vector<int> fN(N,-1);
96  //The order of looping (records or dimensions) does not matter because you have to go through all of them once anyway
97  //Loop over each binned dimension that isn't the first or the last
98  for (unsigned idim=1; idim<=N; idim++) {
99  int nBoundaryPassed = 0;
100  if (fN[idim-1]==-1) fN[idim-1] = mBinBoundaries[idim].size();
101  //Loop over the mRecords vector
102  for (unsigned iRecord=0; iRecord<nRec; iRecord++) {
103  if (mRecords[iRecord].xMin(idim)!=mBinBoundaries[idim][nBoundaryPassed%fN[idim-1]]) {
104  throw cms::Exception("MissingRecord") << "found a missing record in binned dimension " << idim << " after record " << iRecord << std::endl
105  << "\tthe bin lower bound should have been " << mBinBoundaries[idim][nBoundaryPassed%fN[idim-1]]
106  << ", but was instead " << mRecords[iRecord].xMin(idim) << std::endl
107  << "\tall binned dimensions, besides the last one, must have uniform binning." << std::endl
108  << mRecords[iRecord-1] << std::endl << mRecords[iRecord] << std::endl;
109  }
110  else if((iRecord==nRec-1 || mRecords[iRecord].xMin(idim)!=mRecords[iRecord+1].xMin(idim))) {
111  nBoundaryPassed++;
112  }
113  else {
114  continue;
115  }
116  }
117  }
118 }
119 
120 //------------------------------------------------------------------------
121 //--- JetCorrectorParameters::JetCorrectorParametersHelper sanity checks -
122 //--- checks that some conditions are met before finding the bin index ---
123 //------------------------------------------------------------------------
124 void JetCorrectorParametersHelper::binIndexChecks(unsigned N, const std::vector<float>& fX) const
125 {
126  if (N != fX.size())
127  {
128  std::stringstream sserr;
129  sserr<<"The number of binned variables, "<<N<<", doesn't correspond to the number requested, "<<fX.size();
130  handleError("JetCorrectorParametersHelper",sserr.str());
131  }
132 }
133 bool JetCorrectorParametersHelper::binBoundChecks(unsigned dim, const float& value, const float& min, const float& max) const
134 {
135  if (value < min || value > max) return false;
136  else return true;
137 }
138 //------------------------------------------------------------------------
139 //--- JetCorrectorParameters::JetCorrectorParametersHelper binIndexN -----
140 //--- returns the index of the record defined by fX (non-linear search) --
141 //------------------------------------------------------------------------
142 int JetCorrectorParametersHelper::binIndexN(const std::vector<float>& fX,
143  const std::vector<JetCorrectorParameters::Record>& mRecords) const
144 {
145  unsigned Nm1 = SIZE-1;
146  binIndexChecks(SIZE,fX);
147 
148  //Create a container for the indices
149  std::vector<float> fN(SIZE,-1);
150  std::vector<float>::const_iterator tmpIt;
151 
152  // make sure that fX are within the first and last boundaries of mBinBoundaries (other than last dimension)
153  for (unsigned idim=0; idim==0||idim<fX.size()-1; idim++)
154  {
155  if(!binBoundChecks(idim,fX[idim],*mBinBoundaries[idim].begin(),mRecords[size()-1].xMax(idim))) return -1;
156  tmpIt = std::lower_bound(mBinBoundaries[idim].begin(),mBinBoundaries[idim].end(),fX[idim]);
157  // lower_bound finds the entry with the next highest value to fX[0]
158  // so unless the two values are equal, you want the next lowest bin boundary
159  if (tmpIt == mBinBoundaries[idim].end())
160  tmpIt = mBinBoundaries[idim].begin()+mBinBoundaries[idim].size()-1;
161  else if (*tmpIt != fX[idim])
162  tmpIt -= 1;
163  fN[idim] = *tmpIt;
164  }
165 
166  //find the index bounds for the possible values of the last dimension
167  std::pair<size_t,size_t> indexBounds;
168  if(SIZE>1)
169  {
170  tuple_type_Nm1 to_find_Nm1 = gen_tuple<JetCorrectorParameters::MAX_SIZE_DIMENSIONALITY-1>([&](size_t i){return (i<Nm1)?fN[i]:-9999;});
171  if (mMap.find(to_find_Nm1)!=mMap.end())
172  indexBounds = mMap.at(to_find_Nm1);
173  else
174  {
175  std::stringstream sserr;
176  sserr<<"couldn't find the index boundaries for dimension "<<Nm1<<std::endl
177  <<"looking for last bin with N-1 values of "<<to_find_Nm1<<std::endl;
178  handleError("JetCorrectorParametersHelper",sserr.str());
179  return -1;
180  }
181 
182  //Check that the requested value is within the bin boundaries for the last dimension
183  if(!binBoundChecks(Nm1,fX[Nm1],mRecords[indexBounds.first].xMin(Nm1),mRecords[indexBounds.second].xMax(Nm1))) return -1;
184  tmpIt = std::lower_bound(mBinBoundaries[Nm1].begin()+indexBounds.first,mBinBoundaries[Nm1].begin()+indexBounds.second,fX[Nm1]);
185  if (*tmpIt != fX[Nm1] && fX[Nm1]<*(mBinBoundaries[Nm1].begin()+indexBounds.second))
186  tmpIt-=1;
187  fN[Nm1] = *tmpIt;
188  }
189 
190  tuple_type to_find = gen_tuple<JetCorrectorParameters::MAX_SIZE_DIMENSIONALITY>([&](size_t i){return (i<SIZE)?fN[i]:-9999;});
191  return (mIndexMap.find(to_find)!=mIndexMap.end()) ? mIndexMap.at(to_find) : -1;
192 }
Definition: start.py:1
int binIndexN(const std::vector< float > &fX, const std::vector< JetCorrectorParameters::Record > &mRecords) const
void binIndexChecks(unsigned N, const std::vector< float > &fX) const
void checkMiddleBinUniformity(const std::vector< JetCorrectorParameters::Record > &mRecords) const
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
Definition: FindCaloHit.cc:20
std::unordered_map< tuple_type, size_t > mIndexMap
typename generate_tuple_type< float, JetCorrectorParameters::MAX_SIZE_DIMENSIONALITY >::type tuple_type
static const int MAX_SIZE_DIMENSIONALITY
#define end
Definition: vmac.h:39
Definition: value.py:1
T min(T a, T b)
Definition: MathUtil.h:58
typename generate_tuple_type< float, JetCorrectorParameters::MAX_SIZE_DIMENSIONALITY-1 >::type tuple_type_Nm1
int k[5][pyjets_maxn]
bool binBoundChecks(unsigned dim, const float &value, const float &min, const float &max) const
#define N
Definition: blowfish.cc:9
void init(const JetCorrectorParameters::Definitions &mDefinitions, const std::vector< JetCorrectorParameters::Record > &mRecords)
#define begin
Definition: vmac.h:32
std::vector< std::vector< float > > mBinBoundaries
std::unordered_map< tuple_type_Nm1, std::pair< size_t, size_t > > mMap