TauJet/src/TauJetCorrector.cc

0001 #include "JetMETCorrections/TauJet/interface/TauJetCorrector.h"
0002
0003 #include "FWCore/ParameterSet/interface/FileInPath.h"
0004 #include "FWCore/Framework/interface/Event.h"
0005 #include "DataFormats/Common/interface/Handle.h"
0006 #include "FWCore/Framework/interface/EventSetup.h"
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009
0010 #include <vector>
0011 #include <fstream>
0012 #include <sstream>
0013 using namespace std;
0014
0015 double TauJetCorrector::ParametrizationTauJet::value(double et, double eta) const {
0016   double x = et;
0017   double etnew(et);
0018   double etabound = (*theEtabound.find(type)).second;
0019   std::vector<double> taus = (*theParam.find(type)).second;
0020
0021   if (fabs(eta) > etabound) {
0022     //cout << " ===> ETA outside of range - CORRECTION NOT DONE ***" << endl;
0023     //cout << " eta = " << eta <<" pt = " << et << " etabound "<<etabound<<endl;
0024     return et;
0025   }
0026   //cout << "correction parameters: " << taus[0] << " " << taus[1] << " " << taus[2] << endl;
0027   switch (type) {
0028     case 1: {
0029       etnew = 2 * (x - taus[0]) / (taus[1] + sqrt(taus[1] * taus[1] - 4 * taus[0] * taus[2] + 4 * x * taus[2]));
0030       break;
0031     }
0032     case 2: {
0033       etnew = 2 * (x - taus[0]) / (taus[1] + sqrt(taus[1] * taus[1] - 4 * taus[0] * taus[2] + 4 * x * taus[2]));
0034       break;
0035     }
0036     case 3: {
0037       etnew = 2 * (x - taus[0]) / (taus[1] + sqrt(taus[1] * taus[1] - 4 * taus[0] * taus[2] + 4 * x * taus[2]));
0038       break;
0039     }
0040
0041     default:
0042       edm::LogError("TauJetCorrector: Error: unknown parametrization type ")
0043           << type << " in TauJetCorrector. No correction applied" << endl;
0044       //cerr<<"TauJetCorrector: Error: unknown parametrization type '"<<type<<"' in TauJetCorrector. No correction applied"<<endl;
0045       break;
0046   }
0047   return etnew;
0048 }
0049
0050 class JetCalibrationParameterSetTauJet {
0051 public:
0052   JetCalibrationParameterSetTauJet(string tag);
0053   int neta() { return etavector.size(); }
0054   double eta(int ieta) { return etavector[ieta]; }
0055   int type(int ieta) { return typevector[ieta]; }
0056   const vector<double>& parameters(int ieta) { return pars[ieta]; }
0057   bool valid() { return !etavector.empty(); }
0058
0059 private:
0060   vector<double> etavector;
0061   vector<int> typevector;
0062   vector<vector<double> > pars;
0063 };
0064 JetCalibrationParameterSetTauJet::JetCalibrationParameterSetTauJet(string tag) {
0065   std::string file = "JetMETCorrections/TauJet/data/" + tag + ".txt";
0066
0067   edm::FileInPath f1(file);
0068
0069   std::ifstream in((f1.fullPath()).c_str());
0070
0071   //  if ( f1.isLocal() ){
0072   //cout << " Start to read file "<<file<<endl;
0073   string line;
0074   while (std::getline(in, line)) {
0075     if (line.empty() || line[0] == '#')
0076       continue;
0077     istringstream linestream(line);
0078     double par;
0079     int type;
0080     linestream >> par >> type;
0081
0082     //cout<<" Parameter eta = "<<par<<" Type= "<<type<<endl;
0083
0084     etavector.push_back(par);
0085     typevector.push_back(type);
0086     pars.push_back(vector<double>());
0087     while (linestream >> par)
0088       pars.back().push_back(par);
0089   }
0090   //  }
0091   //  else
0092   //    cout<<"The file \""<<file<<"\" was not found in path \""<<f1.fullPath()<<"\"."<<endl;
0093 }
0094
0095 TauJetCorrector::TauJetCorrector(const edm::ParameterSet& fConfig) {
0096   type = fConfig.getParameter<int>("TauTriggerType");
0097   setParameters(fConfig.getParameter<std::string>("tagName"), type);
0098 }
0099
0100 TauJetCorrector::~TauJetCorrector() {
0101   for (ParametersMap::iterator ip = parametrization.begin(); ip != parametrization.end(); ip++)
0102     delete ip->second;
0103 }
0104
0105 void TauJetCorrector::setParameters(std::string aCalibrationType, int itype) {
0106   //cout<< " Start to set parameters "<<endl;
0107   type = itype;
0108   JetCalibrationParameterSetTauJet pset(aCalibrationType);
0109
0110   if ((!pset.valid()) && (aCalibrationType != "no")) {
0111     edm::LogError("TauJetCorrector:Jet Corrections not found ")
0112         << aCalibrationType
0113         << " not found! Cannot apply any correction ... For JetPlusTrack calibration only radii 0.5 and 0.7 are "
0114            "included for JetParton"
0115         << endl;
0116     return;
0117   }
0118   if (aCalibrationType == "no")
0119     return;
0120
0121   map<int, vector<double> > pq;
0122   map<int, vector<double> > pg;
0123   map<int, vector<double> > pqcd;
0124   map<int, double> etaboundx;
0125   int iq = 0;
0126   int ig = 0;
0127   int iqcd = 0;
0128   int mtype = 0;
0129   for (int ieta = 0; ieta < pset.neta(); ieta++) {
0130     if (pset.type(ieta) == 1) {
0131       pq[iq] = pset.parameters(ieta);
0132       iq++;
0133       mtype = (int)(pset.type(ieta));
0134     }
0135     if (pset.type(ieta) == 2) {
0136       pg[ig] = pset.parameters(ieta);
0137       ig++;
0138       mtype = (int)(pset.type(ieta));
0139     }
0140     if (pset.type(ieta) == 3) {
0141       pqcd[iqcd] = pset.parameters(ieta);
0142       iqcd++;
0143       mtype = (int)(pset.type(ieta));
0144     }
0145     if (pset.type(ieta) == -1) {
0146       etaboundx[mtype - 1] = pset.eta(ieta);
0147     }
0148   }
0149
0150   //cout<<" Number of parameters "<<iq<<" "<<ig<<" "<<iqcd<<endl;
0151   int mynum = 0;
0152   for (int ieta = 0; ieta < pset.neta(); ieta++) {
0153     //cout<<" New parmetrization "<<ieta<<" "<<pset.type(ieta)<<endl;
0154
0155     if (pset.type(ieta) == -1)
0156       continue;
0157     if (ieta < iq + 1) {
0158       parametrization[pset.eta(ieta)] =
0159           new ParametrizationTauJet(pset.type(ieta), (*pq.find(ieta)).second, (*etaboundx.find(0)).second);
0160       //cout<<" ALL "<<ieta<<" "<<((*pq.find(ieta)).second)[0]<<" "<<((*pq.find(ieta)).second)[1]<<" "<<
0161       //((*pq.find(ieta)).second)[2]<<endl;
0162     }
0163     if (ieta > iq && ieta < iq + ig + 2) {
0164       mynum = ieta - iq - 1;
0165       parametrization[pset.eta(ieta)] =
0166           new ParametrizationTauJet(pset.type(ieta), (*pg.find(mynum)).second, (*etaboundx.find(1)).second);
0167       //cout<<" One prong "<<((*pg.find(mynum)).second)[0]<<" "<<((*pg.find(mynum)).second)[1]<<" "<<
0168       //((*pg.find(mynum)).second)[2]<<endl;
0169     }
0170     if (ieta > iq + ig + 1) {
0171       mynum = ieta - iq - ig - 2;
0172       //cout<<" Mynum "<<mynum<<" "<<ieta<<" "<<pset.type(ieta)<<endl;
0173       parametrization[pset.eta(ieta)] =
0174           new ParametrizationTauJet(pset.type(ieta), (*pqcd.find(mynum)).second, (*etaboundx.find(2)).second);
0175       //cout<<" Two prongs "<<((*pqcd.find(mynum)).second)[0]<<" "<<((*pqcd.find(mynum)).second)[1]<<" "<<
0176       //((*pqcd.find(mynum)).second)[2]<<endl;
0177     }
0178   }
0179   //cout<<" Parameters inserted into mAlgorithm "<<endl;
0180 }
0181
0182 double TauJetCorrector::correction(const LorentzVector& fJet) const {
0183   //cout<<" Start Apply Corrections "<<endl;
0184   if (parametrization.empty()) {
0185     return 1.;
0186   }
0187
0188   double et = fJet.Et();
0189   double eta = fabs(fJet.Eta());
0190
0191   //cout<<" Et and eta of jet "<<et<<" "<<eta<<endl;
0192
0193   double etnew;
0194   std::map<double, ParametrizationTauJet*>::const_iterator ip = parametrization.upper_bound(eta);
0195   etnew = (--ip)->second->value(et, eta);
0196
0197   //cout<<" The new energy found "<<etnew<<" "<<et<<endl;
0198
0199   float mScale = etnew / et;
0200
0201   return mScale;
0202 }
0203
0204 double TauJetCorrector::correction(const reco::Jet& fJet) const { return correction(fJet.p4()); }