JetAlgorithms/src/QGLikelihoodCalculator.cc

0001 #include "RecoJets/JetAlgorithms/interface/QGLikelihoodCalculator.h"
0002 #include "CondFormats/JetMETObjects/interface/QGLikelihoodObject.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include <cmath>
0005
0006 /// Compute likelihood for a jet using the QGLikelihoodObject information and a set of variables
0007 float QGLikelihoodCalculator::computeQGLikelihood(
0008     const QGLikelihoodObject &QGLParamsColl, float pt, float eta, float rho, std::vector<float> vars) const {
0009   if (!isValidRange(pt, rho, eta, QGLParamsColl.qgValidRange))
0010     return -1;
0011
0012   float Q = 1., G = 1.;
0013   for (unsigned int varIndex = 0; varIndex < vars.size(); ++varIndex) {
0014     auto quarkEntry = findEntry(QGLParamsColl.data, eta, pt, rho, 0, varIndex);
0015     auto gluonEntry = findEntry(QGLParamsColl.data, eta, pt, rho, 1, varIndex);
0016     if (!quarkEntry || !gluonEntry)
0017       return -2;
0018
0019     int binQ = quarkEntry->histogram.findBin(vars[varIndex]);
0020     float Qi = quarkEntry->histogram.binContent(binQ);
0021
0022     int binG = gluonEntry->histogram.findBin(vars[varIndex]);
0023     float Gi = gluonEntry->histogram.binContent(binG);
0024
0025     Q *= Qi;
0026     G *= Gi;
0027   }
0028
0029   if (Q == 0)
0030     return 0;
0031   return Q / (Q + G);
0032 }
0033
0034 /// Find matching entry in vector for a given eta, pt, rho, qgIndex and varIndex
0035 const QGLikelihoodObject::Entry *QGLikelihoodCalculator::findEntry(std::vector<QGLikelihoodObject::Entry> const &data,
0036                                                                    float eta,
0037                                                                    float pt,
0038                                                                    float rho,
0039                                                                    int qgIndex,
0040                                                                    int varIndex) const {
0041   QGLikelihoodParameters myParameters;
0042   myParameters.Rho = rho;
0043   myParameters.Pt = pt;
0044   myParameters.Eta = fabs(eta);
0045   myParameters.QGIndex = qgIndex;
0046   myParameters.VarIndex = varIndex;
0047
0048   auto myDataObject = data.begin();
0049   while (!(myParameters == myDataObject->category)) {
0050     ++myDataObject;
0051     if (myDataObject == data.end()) {
0052       edm::LogWarning("QGLCategoryNotFound")
0053           << "Jet passed qgValidRange criteria, but no category found with rho=" << rho << ", pt=" << pt
0054           << ", eta=" << eta << "\nPlease contact cms-qg-workinggroup@cern.ch" << std::endl;
0055       return nullptr;
0056     }
0057   }
0058   return &*myDataObject;
0059 }
0060
0061 /// Check the valid range of this qg tagger
0062 bool QGLikelihoodCalculator::isValidRange(float pt,
0063                                           float rho,
0064                                           float eta,
0065                                           const QGLikelihoodCategory &qgValidRange) const {
0066   if (pt < qgValidRange.PtMin)
0067     return false;
0068   if (pt > qgValidRange.PtMax)
0069     return false;
0070   if (rho < qgValidRange.RhoMin)
0071     return false;
0072   if (rho > qgValidRange.RhoMax)
0073     return false;
0074   if (fabs(eta) < qgValidRange.EtaMin)
0075     return false;
0076   if (fabs(eta) > qgValidRange.EtaMax)
0077     return false;
0078   return true;
0079 }
0080
0081 /// Return the smeared qgLikelihood value, given input x0 and parameters a, b, min and max
0082 float QGLikelihoodCalculator::smearingFunction(float x0, float a, float b, float min, float max) const {
0083   float x = (x0 - min) / (max - min);
0084   if (x < 0.)
0085     x = 0.;
0086   if (x > 1.)
0087     x = 1.;
0088
0089   float x1 = tanh(a * atanh(2. * x - 1.) + b) / 2. + .5;
0090   if (x <= 0.)
0091     x1 = 0.;
0092   if (x >= 1.)
0093     x1 = 1.;
0094
0095   return x1 * (max - min) + min;
0096 }
0097
0098 // Get systematic smearing
0099 float QGLikelihoodCalculator::systematicSmearing(const QGLikelihoodSystematicsObject &QGLSystematicsColl,
0100                                                  float pt,
0101                                                  float eta,
0102                                                  float rho,
0103                                                  float qgValue,
0104                                                  int qgIndex) const {
0105   if (qgValue < 0 || qgValue > 1)
0106     return -1.;
0107
0108   QGLikelihoodParameters myParameters;
0109   myParameters.Rho = rho;
0110   myParameters.Pt = pt;
0111   myParameters.Eta = fabs(eta);
0112   myParameters.QGIndex = qgIndex;
0113   myParameters.VarIndex = -1;
0114
0115   auto myDataObject = QGLSystematicsColl.data.begin();
0116   while (!(myParameters == myDataObject->systCategory)) {
0117     ++myDataObject;
0118     if (myDataObject == QGLSystematicsColl.data.end())
0119       return -1;  //Smearing not available in the whole qgValidRange: do not throw warnings or errors
0120   }
0121   return smearingFunction(qgValue, myDataObject->a, myDataObject->b, myDataObject->lmin, myDataObject->lmax);
0122 }