CMS 3D CMS Logo

L1GctJetFinderParams.cc
Go to the documentation of this file.
2 
3 #include <iostream>
4 #include <iomanip>
5 #include <math.h>
6 
8 
11 
12 using std::ios;
13 
14 const unsigned L1GctJetFinderParams::NUMBER_ETA_VALUES = 11;
16 
18  rgnEtLsb_(0.),
19  htLsb_(0.),
20  cenJetEtSeed_(0.),
21  forJetEtSeed_(0.),
22  tauJetEtSeed_(0.),
23  tauIsoEtThreshold_(0.),
24  htJetEtThreshold_(0.),
25  mhtJetEtThreshold_(0.),
26  cenForJetEtaBoundary_(0),
27  corrType_(0),
28  jetCorrCoeffs_(),
29  tauCorrCoeffs_(),
30  convertToEnergy_(false),
31  energyConversionCoeffs_()
32 { }
33 
35  double htLsb,
36  double cJetSeed,
37  double fJetSeed,
38  double tJetSeed,
39  double tauIsoEtThresh,
40  double htJetEtThresh,
41  double mhtJetEtThresh,
42  unsigned etaBoundary,
43  unsigned corrType,
44  const std::vector< std::vector<double> >& jetCorrCoeffs,
45  const std::vector< std::vector<double> >& tauCorrCoeffs,
46  bool convertToEnergy,
47  const std::vector<double>& energyConvCoeffs) :
48  rgnEtLsb_(rgnEtLsb),
49  htLsb_(htLsb),
50  cenJetEtSeed_(cJetSeed),
51  forJetEtSeed_(fJetSeed),
52  tauJetEtSeed_(tJetSeed),
53  tauIsoEtThreshold_(tauIsoEtThresh),
54  htJetEtThreshold_(htJetEtThresh),
55  mhtJetEtThreshold_(mhtJetEtThresh),
56  cenForJetEtaBoundary_(etaBoundary),
57  corrType_(corrType),
58  jetCorrCoeffs_(jetCorrCoeffs),
59  tauCorrCoeffs_(tauCorrCoeffs),
60  convertToEnergy_(convertToEnergy),
61  energyConversionCoeffs_(energyConvCoeffs)
62 {
63  // check number of eta bins
64  if (jetCorrCoeffs_.size() != NUMBER_ETA_VALUES ||
67 
68  LogDebug("L1-O2O") << "GCT jet corrections constructed with " << jetCorrCoeffs_.size() << " bins, expected " << NUMBER_ETA_VALUES << std::endl;
69  LogDebug("L1-O2O") << "GCT tau corrections constructed with " << tauCorrCoeffs_.size() << " bins, expected " << N_CENTRAL_ETA_VALUES << std::endl;
70  LogDebug("L1-O2O") << "GCT energy corrections constructed with " << energyConversionCoeffs_.size() << " bins, expected " << NUMBER_ETA_VALUES << std::endl;
71 
72  throw cms::Exception("InconsistentConfig") << "L1GctJetFinderParams constructed with wrong number of eta bins : " << jetCorrCoeffs_.size() << " jets, " << tauCorrCoeffs_.size() << " taus, " << energyConversionCoeffs_.size() << " energy conversion bins" << std::endl;
73 
74  }
75 
76  // check number of coefficients against expectation
77  unsigned expCoeffs = 0;
78  if (corrType_ == 2) expCoeffs=8; // ORCA style
79  if (corrType_ == 3) expCoeffs=4; // Simple
80  if (corrType_ == 4) expCoeffs=15; // piecewise-cubic
81  if (corrType_ == 5) expCoeffs=6; // PF
82 
83  // correction types 1 and 4 can have any number of parameters
84  if (corrType_ != 1 &&
85  corrType_ != 4) {
86  std::vector< std::vector<double> >::const_iterator itr;
87  for (itr=jetCorrCoeffs_.begin(); itr!=jetCorrCoeffs_.end(); ++itr) {
88  if (itr->size() != expCoeffs) {
89  throw cms::Exception("InconsistentConfig") << "L1GctJetFinderParams constructed with " << itr->size() << " jet correction coefficients, when " << expCoeffs << " expected" << std::endl;
90  }
91  }
92  for (itr=tauCorrCoeffs_.begin(); itr!=tauCorrCoeffs_.end(); ++itr) {
93  if (itr->size() != expCoeffs) {
94  throw cms::Exception("InconsistentConfig") << "L1GctJetFinderParams constructed with " << itr->size() << " tau correction coefficients, when " << expCoeffs << " expected"<< std::endl;
95  }
96  }
97  }
98 
99 }
100 
101 
103 
104 //---------------------------------------------------------------------------------------------
105 //
106 // set methods
107 //
108 
109 void L1GctJetFinderParams::setRegionEtLsb (const double rgnEtLsb)
110 {
111  rgnEtLsb_ = rgnEtLsb;
112 }
113 
115  const double fJetSeed,
116  const double tJetSeed,
117  const unsigned etaBoundary)
118 {
119  cenJetEtSeed_ = cJetSeed;
120  forJetEtSeed_ = fJetSeed;
121  tauJetEtSeed_ = tJetSeed;
122  cenForJetEtaBoundary_ = etaBoundary;
123 }
124 
126  const std::vector< std::vector<double> >& jetCorrCoeffs,
127  const std::vector< std::vector<double> >& tauCorrCoeffs)
128 {
129  corrType_ = corrType;
130  jetCorrCoeffs_ = jetCorrCoeffs;
131  tauCorrCoeffs_ = tauCorrCoeffs;
132 }
133 
134 void L1GctJetFinderParams::setJetEtConvertToEnergyOn(const std::vector<double>& energyConvCoeffs)
135 {
136  convertToEnergy_ = true;
137  energyConversionCoeffs_ = energyConvCoeffs;
138 }
139 
141 {
142  convertToEnergy_ = false;
143  energyConversionCoeffs_.clear();
144 }
145 
146 void L1GctJetFinderParams::setHtSumParams(const double htLsb,
147  const double htJetEtThresh,
148  const double mhtJetEtThresh)
149 {
150  htLsb_ = htLsb;
151  htJetEtThreshold_ = htJetEtThresh;
152  mhtJetEtThreshold_ = mhtJetEtThresh;
153 }
154 
155 void L1GctJetFinderParams::setTauAlgorithmParams(const double tauIsoEtThresh)
156 {
157  tauIsoEtThreshold_ = tauIsoEtThresh;
158 }
159 
160 void L1GctJetFinderParams::setParams(const double rgnEtLsb,
161  const double htLsb,
162  const double cJetSeed,
163  const double fJetSeed,
164  const double tJetSeed,
165  const double tauIsoEtThresh,
166  const double htJetEtThresh,
167  const double mhtJetEtThresh,
168  const unsigned etaBoundary,
169  const unsigned corrType,
170  const std::vector< std::vector<double> >& jetCorrCoeffs,
171  const std::vector< std::vector<double> >& tauCorrCoeffs)
172 {
173  setRegionEtLsb (rgnEtLsb);
174  setSlidingWindowParams(cJetSeed, fJetSeed, tJetSeed, etaBoundary);
175  setJetEtCalibrationParams(corrType, jetCorrCoeffs, tauCorrCoeffs);
176  setHtSumParams(htLsb, htJetEtThresh, mhtJetEtThresh);
177  setTauAlgorithmParams(tauIsoEtThresh);
178 }
179 
180 //---------------------------------------------------------------------------------------------
181 //
182 // Jet Et correction methods
183 //
184 
186  const unsigned eta,
187  const bool tauVeto) const {
188 
189  if (eta>=NUMBER_ETA_VALUES) {
190  return 0;
191  } else {
192  double result=0;
193  if ((eta>=cenForJetEtaBoundary_) || tauVeto) {
194  // Use jetCorrCoeffs for central and forward jets.
195  // In forward eta bins we ignore the tau flag (as in the firmware)
196  result=correctionFunction(et, jetCorrCoeffs_.at(eta));
197  } else {
198  // Use tauCorrCoeffs for tau jets (in central eta bins)
199  result=correctionFunction(et, tauCorrCoeffs_.at(eta));
200  }
201  if (convertToEnergy_) { result *= energyConversionCoeffs_.at(eta); }
202  return result;
203  }
204 
205 }
206 
207 
209 uint16_t L1GctJetFinderParams::correctedEtGct(const double correctedEt) const
210 {
211  double scaledEt = correctedEt / htLsb_;
212 
213  uint16_t jetEtOut = static_cast<uint16_t>(scaledEt);
214 
217  } else {
218  return jetEtOut;
219  }
220 }
221 
222 
223 
224 // private methods
225 
226 double L1GctJetFinderParams::correctionFunction(const double Et, const std::vector<double>& coeffs) const
227 {
228  double result=0;
229  switch (corrType_)
230  {
231  case 0: // no correction
232  result = Et;
233  break;
234  case 1: // power series correction
235  result = powerSeriesCorrect(Et, coeffs);
236  break;
237  case 2: // ORCA style correction
238  result = orcaStyleCorrect(Et, coeffs);
239  break;
240  case 3: // simple correction (JetMEt style)
241  result = simpleCorrect(Et, coeffs);
242  break;
243  case 4: // piecwise cubic correction a la Greg Landsberg et al
244  result = piecewiseCubicCorrect(Et, coeffs);
245  break;
246  case 5: // PF correction
247  result = pfCorrect(Et, coeffs);
248  break;
249  default:
250  result = Et;
251  }
252  return result;
253 }
254 
255 double L1GctJetFinderParams::powerSeriesCorrect(const double Et, const std::vector<double>& coeffs) const
256 {
257  double corrEt = Et;
258  for (unsigned i=0; i<coeffs.size();i++) {
259  corrEt += coeffs.at(i)*pow(Et,(int)i);
260  }
261  return corrEt;
262 }
263 
264 double L1GctJetFinderParams::orcaStyleCorrect(const double Et, const std::vector<double>& coeffs) const
265 {
266 
267  // The coefficients are arranged in groups of four
268  // The first in each group is a threshold value of Et.
269 
270  std::vector<double>::const_iterator next_coeff=coeffs.begin();
271  while (next_coeff != coeffs.end()) {
272  double threshold = *next_coeff++;
273  double A = *next_coeff++;
274  double B = *next_coeff++;
275  double C = *next_coeff++;
276  if (Et>threshold) {
277  // This function is an inverse quadratic:
278  // (input Et) = A + B*(output Et) + C*(output Et)^2
279  return 2*(Et-A)/(B+sqrt(B*B-4*A*C+4*Et*C));
280  }
281  // If we are below all specified thresholds (or the vector is empty), return output=input.
282  }
283  return Et;
284 }
285 
286 double L1GctJetFinderParams::simpleCorrect(const double Et, const std::vector<double>& coeffs) const
287 {
288  // function is :
289  // Et_out = A + B/[ (log10(Et_in))^C + D ] + E/Et_in
290  //
291  // fitcor_as_str = "[0]+[1]/(pow(log10(x),[2])+[3]) + [4]/x"
292  //
293 
294  return coeffs.at(0) + coeffs.at(1)/(pow(log10(Et),coeffs.at(2))+coeffs.at(3)) + (coeffs.at(4)/Et);
295 
296 }
297 
298 double L1GctJetFinderParams::piecewiseCubicCorrect(const double Et, const std::vector<double>& coeffs) const
299 {
300 
301  // The correction fuction is a set of 3rd order polynomials
302  // Et_out = Et_in + (p0 + p1*Et_in + p2*Et_in^2 + p3*Et_in^3)
303  // with different coefficients for different energy ranges.
304  // The parameters are arranged in groups of five.
305  // The first in each group is a threshold value of input Et,
306  // followed by the four coefficients for the cubic function.
307  double etOut = Et;
308  std::vector<double>::const_iterator next_coeff=coeffs.begin();
309  while (next_coeff != coeffs.end()) {
310 
311  // Read the coefficients from the vector
312  double threshold = *next_coeff++;
313  double A = *next_coeff++; //p0
314  double B = *next_coeff++; //p1
315  double C = *next_coeff++; //p2
316  double D = *next_coeff++; //p3
317 
318  // Check we are in the right energy range and make correction
319  if (Et>threshold) {
320  etOut += (A + etOut*(B + etOut*(C + etOut*D))) ;
321  break;
322  }
323 
324  }
325  return etOut;
326 
327 }
328 
329 double L1GctJetFinderParams::pfCorrect(const double et, const std::vector<double>& coeffs) const
330 {
331  //
332  // corr_factor = [0]+[1]/(pow(log10(x),2)+[2])+[3]*exp(-[4]*(log10(x)-[5])*(log10(x)-[5]))
333  // Et_out = Et_in * corr_factor
334  //
335 
336  return et * (coeffs.at(0)+coeffs.at(1)/(pow(log10(et),2)+coeffs.at(2))+coeffs.at(3)*exp(-coeffs.at(4)*(log10(et)-coeffs.at(5))*(log10(et)-coeffs.at(5))));
337 
338 }
339 
340 
341 std::ostream& operator << (std::ostream& os, const L1GctJetFinderParams& fn)
342 {
343  // os << std::setprecision(2);
344 
345  os << "=== Level-1 GCT : Jet Finder Parameters ===" << std::endl;
346  os << "RCT region LSB : " << std::fixed << fn.getRgnEtLsbGeV() << " GeV" << std::endl;
347  os << "Central jet seed threshold : " << std::fixed << fn.getCenJetEtSeedGeV() << " GeV" << std::endl;
348  os << "Tau jet seed threshold : " << std::fixed << fn.getTauJetEtSeedGeV() << " GeV" << std::endl;
349  os << "Forward jet seed threshold : " << std::fixed << fn.getForJetEtSeedGeV() << " GeV" << std::endl;
350  os << "Tau isolation threshold : " << std::fixed << fn.getTauIsoEtThresholdGeV() << " GeV" << std::endl;
351  os << "Ht jet Et threshold : " << std::fixed << fn.getHtJetEtThresholdGeV() << " GeV" << std::endl;
352  os << "MHt jet Et threshold : " << std::fixed << fn.getMHtJetEtThresholdGeV() << " GeV" << std::endl;
353  os << "Ht LSB : " << std::fixed << fn.getHtLsbGeV() << " GeV" << std::endl;
354  os << "Central/Forward boundary : " << std::fixed << fn.getCenForJetEtaBoundary() << std::endl;
355 
356  os << std::endl;
357 
358  os << std::setprecision(6);
359  os << ios::scientific;
360 
361  os << "=== Level-1 GCT : Jet Et Calibration Function ===" << std::endl;
362  if (fn.getCorrType() == 0) {
363  os << "No jet energy corrections applied" << std::endl;
364  } else {
365  switch (fn.getCorrType())
366  {
367  case 1:
368  os << "Function = Power series" << std::endl;
369  break;
370  case 2:
371  os << "Function = ORCA" << std::endl;
372  break;
373  case 3:
374  os << "Function = Simple" << std::endl;
375  break;
376  case 4:
377  os << "Function = PiecewiseCubic" << std::endl;
378  break;
379  case 5:
380  os << "Function = PF" << std::endl;
381  break;
382  default:
383  os << "Unrecognised" << std::endl;
384  break;
385  }
386  std::vector< std::vector<double> > jetCoeffs = fn.getJetCorrCoeffs();
387  std::vector< std::vector<double> > tauCoeffs = fn.getTauCorrCoeffs();
388 
389  os << "Non-tau jet correction coefficients" << std::endl;
390  for (unsigned i=0; i<jetCoeffs.size(); i++){
391  os << "Eta =" << std::setw(2) << i;
392  if (jetCoeffs.at(i).empty()) {
393  os << ", no coefficients";
394  } else {
395  os << " Coefficients = ";
396  for (unsigned j=0; j<jetCoeffs.at(i).size();j++){
397  os << jetCoeffs.at(i).at(j) << ", ";
398  }
399  }
400  os << std::endl;
401  }
402  os << "Tau jet correction coefficients" << std::endl;
403  for (unsigned i=0; i<tauCoeffs.size(); i++){
404  os << "Eta =" << std::setw(2) << i;
405  if (tauCoeffs.at(i).empty()) {
406  os << ", no coefficients";
407  } else {
408  os << " Coefficients = ";
409  for (unsigned j=0; j<tauCoeffs.at(i).size();j++){
410  os << tauCoeffs.at(i).at(j) << ", ";
411  }
412  }
413  os << std::endl;
414  }
415  }
416 
417  os.unsetf(ios::fixed | ios::scientific);
418 
419  return os;
420 }
421 
422 
#define LogDebug(id)
double powerSeriesCorrect(const double Et, const std::vector< double > &coeffs) const
static const unsigned N_CENTRAL_ETA_VALUES
Number of eta bins used in correction.
void setHtSumParams(const double htLsb, const double htJetEtThresh, const double mhtJetEtThresh)
std::vector< std::vector< double > > jetCorrCoeffs_
void setRegionEtLsb(const double rgnEtLsb)
std::vector< std::vector< double > > tauCorrCoeffs_
const std::vector< std::vector< double > > & getJetCorrCoeffs() const
double getHtLsbGeV() const
uint16_t correctedEtGct(const double correctedEt) const
Convert the corrected Et value to a linear Et for Ht summing.
double getCenJetEtSeedGeV() const
double getTauJetEtSeedGeV() const
double correctedEtGeV(const double et, const unsigned eta, const bool tauVeto) const
Eta takes a value from 0-10, corresponding to jet regions running from eta=0.0 to eta=5...
void setParams(const double rgnEtLsb, const double htLsb, const double cJetSeed, const double fJetSeed, const double tJetSeed, const double tauIsoEtThresh, const double htJetEtThresh, const double mhtJetEtThresh, const unsigned etaBoundary, const unsigned corrType, const std::vector< std::vector< double > > &jetCorrCoeffs, const std::vector< std::vector< double > > &tauCorrCoeffs)
double correctionFunction(const double Et, const std::vector< double > &coeffs) const
unsigned getCorrType() const
Access to jet Et calibration parameters.
T sqrt(T t)
Definition: SSEVec.h:18
double simpleCorrect(const double Et, const std::vector< double > &coeffs) const
void setJetEtConvertToEnergyOn(const std::vector< double > &energyConvCoeffs)
double getMHtJetEtThresholdGeV() const
static const std::string B
static const unsigned jetCalibratedEtMax
double getTauIsoEtThresholdGeV() const
double piecewiseCubicCorrect(const double Et, const std::vector< double > &coeffs) const
static const unsigned NUMBER_ETA_VALUES
Number of eta bins used in correction.
DecomposeProduct< arg, typename Div::arg > D
Definition: Factorize.h:150
unsigned getCenForJetEtaBoundary() const
void setJetEtCalibrationParams(const unsigned corrType, const std::vector< std::vector< double > > &jetCorrCoeffs, const std::vector< std::vector< double > > &tauCorrCoeffs)
et
define resolution functions of each parameter
double getHtJetEtThresholdGeV() const
void setTauAlgorithmParams(const double tauIsoEtThresh)
double pfCorrect(const double Et, const std::vector< double > &coeffs) const
double getRgnEtLsbGeV() const
double orcaStyleCorrect(const double Et, const std::vector< double > &coeffs) const
void setSlidingWindowParams(const double cJetSeed, const double fJetSeed, const double tJetSeed, const unsigned etaBoundary)
double getForJetEtSeedGeV() const
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:40
std::vector< double > energyConversionCoeffs_
std::ostream & operator<<(std::ostream &os, const L1GctJetFinderParams &fn)
Overload << operator.
const std::vector< std::vector< double > > & getTauCorrCoeffs() const