JetParton/src/JetPartonCorrector.cc

0001 #include "JetMETCorrections/JetParton/interface/JetPartonCorrector.h"
0002 #include "FWCore/ParameterSet/interface/FileInPath.h"
0003 #include "FWCore/Framework/interface/Event.h"
0004 #include "DataFormats/Common/interface/Handle.h"
0005 #include "FWCore/Framework/interface/EventSetup.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008 #include <vector>
0009 #include <fstream>
0010 #include <sstream>
0011 using namespace std;
0012
0013 namespace JetPartonNamespace {
0014   class UserPartonMixture {
0015   public:
0016     UserPartonMixture() {}
0017     virtual ~UserPartonMixture() {}
0018     virtual double mixt(double et, double eta) {
0019       // The mixture of quark and gluons corresponds to QCD jets.
0020       double f = 0.2 + et / 1000. * (1. + .13 * exp(eta));
0021       return f;
0022     }
0023   };
0024
0025   class ParametrizationJetParton {
0026   public:
0027     ParametrizationJetParton(int thePartonMixture,
0028                              const vector<double>& x,
0029                              const vector<double>& y,
0030                              const vector<double>& z) {
0031       type = thePartonMixture;
0032       pq = x;
0033       pg = y;
0034       pqcd = z;
0035     }
0036
0037     double value(double arg1, double arg2) const;
0038
0039   private:
0040     int type;
0041     std::vector<double> pq;
0042     std::vector<double> pg;
0043     std::vector<double> pqcd;
0044   };
0045
0046   double ParametrizationJetParton::value(double e, double eta) const {
0047     double enew(e);
0048     double x = e;
0049     double y = eta;
0050
0051     if (abs(x - pq[2]) < pq[1] / pq[0] || abs(x - pg[2]) < pg[1] / pg[0] || abs(x - pqcd[2]) < pqcd[1] / pqcd[0]) {
0052       return enew;
0053     }
0054
0055     double kjetq = pq[0] - pq[1] / (x - pq[2]);
0056     double kjetg = pg[0] - pg[1] / (x - pg[2]);
0057     double kjetqcd = pqcd[0] - pqcd[1] / (x - pqcd[2]);
0058
0059     switch (type) {
0060       case 1: {
0061         if (abs(kjetq) > 0.0001)
0062           enew = e / kjetq;
0063         break;
0064       }
0065       case 2: {
0066         if (abs(kjetg) > 0.0001)
0067           enew = e / kjetg;
0068         break;
0069       }
0070       case 3: {
0071         if (abs(kjetqcd) > 0.0001)
0072           enew = e / kjetqcd;
0073         break;
0074       }
0075       case 4: {
0076         cout << "[Jets] JetPartonCorrector: Warning! Calibration to b-quark - does not implemented yet. Light quark "
0077                 "calibration is used instead "
0078              << endl;
0079         if (abs(kjetq) > 0.0001)
0080           enew = e / kjetq;
0081         break;
0082       }
0083       case 100: {
0084         UserPartonMixture upm;
0085         double f = upm.mixt(x, y);
0086         double kjet = (f * kjetq + kjetg) / (f + 1);
0087         if (abs(kjet) > 0.0001)
0088           enew = e / kjet;
0089         break;
0090       }
0091
0092       default:
0093         cerr << "[Jets] JetPartonCorrector: Error! unknown parametrization type = " << type
0094              << " No correction applied ..." << endl;
0095         break;
0096     }
0097     return enew;
0098   }
0099
0100   class JetPartonCalibrationParameterSet {
0101   public:
0102     JetPartonCalibrationParameterSet(string tag);
0103     int neta() { return etavector.size(); }
0104     double eta(int ieta) { return etavector[ieta]; }
0105     int type(int ieta) { return typevector[ieta]; }
0106     const vector<double>& parameters(int ieta) { return pars[ieta]; }
0107     bool valid() { return !etavector.empty(); }
0108
0109   private:
0110     vector<double> etavector;
0111     vector<int> typevector;
0112     vector<vector<double> > pars;
0113   };
0114
0115   JetPartonCalibrationParameterSet::JetPartonCalibrationParameterSet(string tag) {
0116     std::string file = "JetMETCorrections/JetParton/data/" + tag + ".txt";
0117
0118     edm::FileInPath f1(file);
0119
0120     std::ifstream in((f1.fullPath()).c_str());
0121
0122     //  if ( f1.isLocal() ){
0123     string line;
0124     while (std::getline(in, line)) {
0125       if (line.empty() || line[0] == '#')
0126         continue;
0127       istringstream linestream(line);
0128       double par;
0129       int type;
0130       linestream >> par >> type;
0131       etavector.push_back(par);
0132       typevector.push_back(type);
0133       pars.push_back(vector<double>());
0134       while (linestream >> par)
0135         pars.back().push_back(par);
0136     }
0137     //  }
0138     //  else
0139     //    if (tag!="no") { cout<<"The file \""<<file<<"\" was not found in path \""<<f1.fullPath()<<"\"."<<endl; }
0140   }
0141 }  // namespace JetPartonNamespace
0142
0143 JetPartonCorrector::JetPartonCorrector(const edm::ParameterSet& fConfig) {
0144   thePartonMixture = fConfig.getParameter<int>("MixtureType");
0145   theJetFinderRadius = fConfig.getParameter<double>("Radius");
0146   setParameters(fConfig.getParameter<std::string>("tagName"), theJetFinderRadius, thePartonMixture);
0147 }
0148
0149 JetPartonCorrector::~JetPartonCorrector() {
0150   for (ParametersMap::iterator ip = parametrization.begin(); ip != parametrization.end(); ip++)
0151     delete ip->second;
0152 }
0153
0154 void JetPartonCorrector::setParameters(std::string aCalibrationType, double aJetFinderRadius, int aPartonMixture) {
0155   theJetFinderRadius = aJetFinderRadius;
0156   thePartonMixture = aPartonMixture;
0157
0158   JetPartonNamespace::JetPartonCalibrationParameterSet pset(aCalibrationType);
0159
0160   if ((!pset.valid()) && (aCalibrationType != "no")) {
0161     edm::LogError("JetPartonCorrector: Jet Corrections not found ")
0162         << aCalibrationType
0163         << " not found! Cannot apply any correction ... For JetPlusTrack calibration only radii 0.5 and 0.7 are "
0164            "included for JetParton"
0165         << endl;
0166     return;
0167   }
0168   if (aCalibrationType == "no")
0169     return;
0170
0171   map<int, vector<double> > pq;
0172   map<int, vector<double> > pg;
0173   map<int, vector<double> > pqcd;
0174   int iq = 0;
0175   int ig = 0;
0176   int iqcd = 0;
0177   for (int ieta = 0; ieta < pset.neta(); ieta++) {
0178     if (pset.type(ieta) == 1) {
0179       pq[iq] = pset.parameters(ieta);
0180       iq++;
0181     };
0182     if (pset.type(ieta) == 2) {
0183       pg[ig] = pset.parameters(ieta);
0184       ig++;
0185     };
0186     if (pset.type(ieta) == 3) {
0187       pqcd[iqcd] = pset.parameters(ieta);
0188       iqcd++;
0189     };
0190   }
0191
0192   for (int ieta = 0; ieta < iq; ieta++) {
0193     parametrization[pset.eta(ieta)] = new JetPartonNamespace::ParametrizationJetParton(
0194         thePartonMixture, (*pq.find(ieta)).second, (*pg.find(ieta)).second, (*pqcd.find(ieta)).second);
0195   }
0196 }
0197 double JetPartonCorrector::correction(const LorentzVector& fJet) const {
0198   if (parametrization.empty()) {
0199     return 1.;
0200   }
0201
0202   double et = fJet.Et();
0203   double eta = fabs(fJet.Eta());
0204
0205   //if(eta<10) { eta=abs(fJet.getY()); }
0206
0207   double etnew;
0208   std::map<double, JetPartonNamespace::ParametrizationJetParton*>::const_iterator ip = parametrization.upper_bound(eta);
0209   if (ip == parametrization.begin()) {
0210     etnew = ip->second->value(et, eta);
0211   } else if (ip == parametrization.end()) {
0212     etnew = (--ip)->second->value(et, eta);
0213   } else {
0214     double eta2 = ip->first;
0215     double et2 = ip->second->value(et, eta);
0216     ip--;
0217     double eta1 = ip->first;
0218     double et1 = ip->second->value(et, eta);
0219
0220     etnew = (eta2 * et1 - eta1 * et2 + eta * et2 - eta * et1) / (eta2 - eta1);
0221   }
0222   cout << " JetParton::The new energy found " << etnew << " " << et << endl;
0223   float mScale = 1000.;
0224
0225   if (et > 0.001)
0226     mScale = etnew / et;
0227
0228   return mScale;
0229 }