JetPlusTracks/plugins/SimpleZSPJPTJetCorrector.cc

0001 //
0002 // Original Author:  Fedor Ratnikov Dec 27, 2006
0003 //
0004 // ZSP Jet Corrector
0005 //
0006 #include "SimpleZSPJPTJetCorrector.h"
0007
0008 #include <vector>
0009 #include <iostream>
0010 #include <fstream>
0011 #include <sstream>
0012
0013 #include "Math/PtEtaPhiE4D.h"
0014 #include "Math/LorentzVector.h"
0015 typedef ROOT::Math::LorentzVector<ROOT::Math::PtEtaPhiE4D<double> > PtEtaPhiELorentzVectorD;
0016 typedef ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > XYZTLorentzVectorD;
0017
0018 using namespace std;
0019
0020 SimpleZSPJPTJetCorrector::SimpleZSPJPTJetCorrector() { debug_ = false; }
0021
0022 SimpleZSPJPTJetCorrector::SimpleZSPJPTJetCorrector(const std::string& fDataFile) {
0023   debug_ = false;
0024   mParameters = new JetCorrectorParameters(fDataFile, "");
0025
0026   mFunc = new TFormula("function", ((mParameters->definitions()).formula()).c_str());
0027
0028   // Read parameters
0029   if (debug_) {
0030     std::cout << " Size of parameters as read by SimpleJetCorrectorParameters " << mParameters->size() << std::endl;
0031     for (unsigned int i = 0; i < mParameters->size(); i++) {
0032       const std::vector<float> p = mParameters->record(i).parameters();
0033       for (std::vector<float>::const_iterator j = p.begin(); j < p.end(); j++) {
0034         std::cout << " Parameter number " << mParameters->record(i).xMin(0) << " " << mParameters->record(i).xMax(0)
0035                   << " " << (*j) << std::endl;
0036       }
0037     }
0038   }  // debug
0039 }
0040
0041 SimpleZSPJPTJetCorrector::~SimpleZSPJPTJetCorrector() {}
0042
0043 double SimpleZSPJPTJetCorrector::correctionPtEtaPhiE(double fPt, double fEta, double fPhi, double fE) const {
0044   double costhetainv = cosh(fEta);
0045   return correctionEtEtaPhiP(fE / costhetainv, fEta, fPhi, fPt * costhetainv);
0046 }
0047
0048 double SimpleZSPJPTJetCorrector::correctionEtEtaPhiP(double fEt, double fEta, double fPhi, double fP) const {
0049   double et = fEt;
0050   double eta = fabs(fEta);
0051   double factor = 1.;
0052
0053   // Define Eta bin for parametrization
0054   std::vector<float> xx;
0055   xx.push_back(eta);
0056   int band = mParameters->binIndex(xx);
0057
0058   if (band < 0)
0059     return factor;
0060
0061   const std::vector<float> p = mParameters->record(band).parameters();
0062
0063   // Set parameters
0064   for (unsigned int i = 0; i < p.size(); i++) {
0065     mFunc->SetParameter(i, p[i]);
0066   }
0067
0068   if (debug_) {
0069     cout << " Et and eta of jet " << et << " " << eta << " bin " << band << " " << p[1] << " " << p[2] << " " << p[3]
0070          << " " << p[4] << " " << p[5] << endl;
0071   }
0072
0073   //  double koef = 1. - p[2]*exp(p[3]*et)-p[4]*exp(p[5]*et);
0074
0075   double koef = 1. - mFunc->Eval(et);
0076
0077   // If et calojet less then some value - use correction on the boundary
0078
0079   if (et < p[0])
0080     koef = 1. - mFunc->Eval(p[0]);
0081
0082   //
0083   if (koef <= 0.000001) {
0084     if (debug_)
0085       std::cout << "SimpleZSPJPTJetCorrector::Problem with ZSP corrections " << koef << std::endl;
0086     koef = 1.;
0087   }
0088
0089   double etnew = et / koef;
0090
0091   if (debug_)
0092     cout << " The new energy found " << etnew << " " << et << " " << koef << endl;
0093
0094   return etnew / et;
0095 }
0096
0097 double SimpleZSPJPTJetCorrector::correctionPUEtEtaPhiP(double fEt, double fEta, double fPhi, double fP) const {
0098   double et = fEt;
0099   double eta = fabs(fEta);
0100   double factor = 1.;
0101
0102   // Define Eta bin for parametrization
0103   std::vector<float> xx;
0104   xx.push_back(eta);
0105   int band = mParameters->binIndex(xx);
0106
0107   if (band < 0)
0108     return factor;
0109
0110   const std::vector<float> p = mParameters->record(band).parameters();
0111
0112   if (debug_) {
0113     cout << " Et and eta of jet " << et << " " << eta << " bin " << band << std::endl;
0114   }
0115
0116   double koef = (et - p[2]) / et;
0117   double etnew = et / koef;
0118
0119   if (debug_)
0120     cout << " The new energy found " << etnew << " " << et << endl;
0121
0122   return etnew / et;
0123 }