File indexing completed on 2023-03-17 11:10:46
0001 #include <memory>
0002 #include "TH2F.h"
0003 #include "TGraph.h"
0004 #include "TGraphErrors.h"
0005 #include "TRandom.h"
0006 #include "TLorentzVector.h"
0007
0008 #include "FWCore/Framework/interface/Frameworkfwd.h"
0009 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0010
0011 #include "FWCore/Framework/interface/Event.h"
0012 #include "FWCore/Framework/interface/MakerMacros.h"
0013 #include "FWCore/Framework/interface/EventSetup.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/Framework/interface/ESHandle.h"
0016 #include "JetMETCorrections/JetCorrector/interface/JetCorrector.h"
0017 #include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h"
0018 #include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h"
0019 #include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
0020
0021 #include "FWCore/ServiceRegistry/interface/Service.h"
0022 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0023
0024
0025
0026
0027
0028 class JetCorrectorDemo : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0029 public:
0030 explicit JetCorrectorDemo(const edm::ParameterSet &);
0031 ~JetCorrectorDemo() override;
0032 typedef reco::Particle::LorentzVector LorentzVector;
0033
0034 private:
0035 void beginJob() override;
0036 void analyze(const edm::Event &, const edm::EventSetup &) override;
0037 void endJob() override;
0038
0039 edm::EDGetTokenT<reco::JetCorrector> mJetCorrector;
0040 std::string mPayloadName, mUncertaintyTag, mUncertaintyFile;
0041 edm::ESGetToken<JetCorrectorParametersCollection, JetCorrectionsRecord> mPayloadToken;
0042 bool mDebug, mUseCondDB;
0043 int mNHistoPoints, mNGraphPoints;
0044 double mEtaMin, mEtaMax, mPtMin, mPtMax;
0045 std::vector<double> mVEta, mVPt;
0046 double vjec_eta[100][1000], vjec_pt[100][1000], vpt[100][1000], vptcor[100][1000], veta[100][1000];
0047 double vjecUnc_eta[100][1000], vUnc_eta[100][1000], vjecUnc_pt[100][1000], vUnc_pt[100][1000], vex_eta[100][1000],
0048 vex_pt[100][1000];
0049 edm::Service<TFileService> fs;
0050 TH2F *mJECvsEta, *mJECvsPt;
0051 TGraphErrors *mVGraphEta[100], *mVGraphPt[100], *mVGraphCorPt[100];
0052 TGraph *mUncEta[100], *mUncCorPt[100];
0053 TRandom *mRandom;
0054 };
0055
0056
0057
0058
0059 JetCorrectorDemo::JetCorrectorDemo(const edm::ParameterSet &iConfig) {
0060 mJetCorrector = consumes<reco::JetCorrector>(iConfig.getParameter<edm::InputTag>("JetCorrector"));
0061 mPayloadName = iConfig.getParameter<std::string>("PayloadName");
0062 mUncertaintyTag = iConfig.getParameter<std::string>("UncertaintyTag");
0063 mUncertaintyFile = iConfig.getParameter<std::string>("UncertaintyFile");
0064 mNHistoPoints = iConfig.getParameter<int>("NHistoPoints");
0065 mNGraphPoints = iConfig.getParameter<int>("NGraphPoints");
0066 mEtaMin = iConfig.getParameter<double>("EtaMin");
0067 mEtaMax = iConfig.getParameter<double>("EtaMax");
0068 mPtMin = iConfig.getParameter<double>("PtMin");
0069 mPtMax = iConfig.getParameter<double>("PtMax");
0070 mVEta = iConfig.getParameter<std::vector<double> >("VEta");
0071 mVPt = iConfig.getParameter<std::vector<double> >("VPt");
0072 mDebug = iConfig.getUntrackedParameter<bool>("Debug", false);
0073 mUseCondDB = iConfig.getUntrackedParameter<bool>("UseCondDB", false);
0074 if (mUseCondDB) {
0075 mPayloadToken = esConsumes(edm::ESInputTag("", mPayloadName));
0076 }
0077 usesResource(TFileService::kSharedResource);
0078 }
0079
0080 JetCorrectorDemo::~JetCorrectorDemo() {}
0081
0082 void JetCorrectorDemo::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0083 edm::Handle<reco::JetCorrector> corrector;
0084 iEvent.getByToken(mJetCorrector, corrector);
0085 std::unique_ptr<JetCorrectionUncertainty> jecUnc;
0086 if (!mUncertaintyTag.empty()) {
0087 if (mUseCondDB) {
0088 const JetCorrectorParametersCollection &JetCorParColl = iSetup.getData(mPayloadToken);
0089 JetCorrectorParameters const &JetCorPar = JetCorParColl[mUncertaintyTag];
0090 jecUnc = std::make_unique<JetCorrectionUncertainty>(JetCorPar);
0091 std::cout << "Configured Uncertainty from CondDB" << std::endl;
0092 } else {
0093 edm::FileInPath fip("CondFormats/JetMETObjects/data/" + mUncertaintyFile + ".txt");
0094 jecUnc = std::make_unique<JetCorrectionUncertainty>(fip.fullPath());
0095 }
0096 }
0097 double jec, rawPt, corPt, eta, unc;
0098 TLorentzVector P4;
0099 double dEta = (mEtaMax - mEtaMin) / mNGraphPoints;
0100 for (int i = 0; i < mNHistoPoints; i++) {
0101 rawPt = mRandom->Uniform(mPtMin, mPtMax);
0102 eta = mRandom->Uniform(mEtaMin, mEtaMax);
0103 P4.SetPtEtaPhiE(rawPt, eta, 0, 0);
0104 LorentzVector rawP4(P4.Px(), P4.Py(), P4.Pz(), P4.E());
0105 jec = corrector->correction(rawP4);
0106 mJECvsEta->Fill(eta, jec);
0107 mJECvsPt->Fill(rawPt, jec);
0108 }
0109
0110 for (unsigned ieta = 0; ieta < mVEta.size(); ieta++) {
0111 double rPt = pow((3500. / TMath::CosH(mVEta[ieta])) / mPtMin, 1. / mNGraphPoints);
0112 for (int i = 0; i < mNGraphPoints; i++) {
0113 rawPt = mPtMin * pow(rPt, i);
0114 eta = mVEta[ieta];
0115 vpt[ieta][i] = rawPt;
0116 P4.SetPtEtaPhiE(rawPt, eta, 0, 0);
0117 LorentzVector rawP4(P4.Px(), P4.Py(), P4.Pz(), P4.E());
0118 jec = corrector->correction(rawP4);
0119 vjec_eta[ieta][i] = jec;
0120 vptcor[ieta][i] = rawPt * jec;
0121 unc = 0.0;
0122 if (!mUncertaintyTag.empty()) {
0123 jecUnc->setJetEta(eta);
0124 jecUnc->setJetPt(rawPt * jec);
0125 unc = jecUnc->getUncertainty(true);
0126 }
0127 vjecUnc_eta[ieta][i] = unc * jec;
0128 vUnc_eta[ieta][i] = unc;
0129 vex_eta[ieta][i] = 0.0;
0130 if (mDebug)
0131 std::cout << rawPt << " " << eta << " " << jec << " " << rawPt * jec << " " << unc << std::endl;
0132 }
0133 }
0134
0135 for (unsigned ipt = 0; ipt < mVPt.size(); ipt++) {
0136 for (int i = 0; i < mNGraphPoints; i++) {
0137 eta = mEtaMin + i * dEta;
0138 corPt = mVPt[ipt];
0139 veta[ipt][i] = eta;
0140
0141 double e = 1.0;
0142 int nLoop(0);
0143 rawPt = corPt;
0144 while (e > 0.0001 && nLoop < 10) {
0145 P4.SetPtEtaPhiE(rawPt, eta, 0, 0);
0146 LorentzVector rawP4(P4.Px(), P4.Py(), P4.Pz(), P4.E());
0147 jec = corrector->correction(rawP4);
0148 double tmp = rawPt * jec;
0149 e = fabs(tmp - corPt) / corPt;
0150 if (jec > 0)
0151 rawPt = corPt / jec;
0152 nLoop++;
0153 }
0154
0155 P4.SetPtEtaPhiE(rawPt, eta, 0, 0);
0156 LorentzVector rawP4(P4.Px(), P4.Py(), P4.Pz(), P4.E());
0157 jec = corrector->correction(rawP4);
0158 vjec_pt[ipt][i] = jec;
0159 unc = 0.0;
0160 if (!mUncertaintyTag.empty()) {
0161 jecUnc->setJetEta(eta);
0162 jecUnc->setJetPt(corPt);
0163 unc = jecUnc->getUncertainty(true);
0164 }
0165 vjecUnc_pt[ipt][i] = unc * jec;
0166 vUnc_pt[ipt][i] = unc;
0167 vex_pt[ipt][i] = 0.0;
0168 if (mDebug)
0169 std::cout << rawPt << " " << eta << " " << jec << " " << rawPt * jec << " " << unc << std::endl;
0170 }
0171 }
0172 }
0173
0174 void JetCorrectorDemo::beginJob() {
0175 if (mNGraphPoints > 1000)
0176 throw cms::Exception("JetCorrectorDemo", "Too many graph points !!! Maximum is 1000 !!!");
0177 if (mVEta.size() > 100)
0178 throw cms::Exception("JetCorrectorDemo", "Too many eta values !!! Maximum is 100 !!!");
0179 if (mVPt.size() > 100)
0180 throw cms::Exception("JetCorrectorDemo", "Too many pt values !!! Maximum is 100 !!!");
0181 mJECvsEta = fs->make<TH2F>("JECvsEta", "JECvsEta", 200, mEtaMin, mEtaMax, 100, 0, 5);
0182 mJECvsPt = fs->make<TH2F>("JECvsPt", "JECvsPt", 200, mPtMin, mPtMax, 100, 0, 5);
0183 mRandom = new TRandom();
0184 mRandom->SetSeed(0);
0185 }
0186
0187 void JetCorrectorDemo::endJob() {
0188 char name[1000];
0189 for (unsigned ipt = 0; ipt < mVPt.size(); ipt++) {
0190 mVGraphEta[ipt] = fs->make<TGraphErrors>(mNGraphPoints, veta[ipt], vjec_pt[ipt], vex_pt[ipt], vjecUnc_pt[ipt]);
0191 sprintf(name, "JEC_vs_Eta_CorPt%1.1f", mVPt[ipt]);
0192 mVGraphEta[ipt]->SetName(name);
0193 mUncEta[ipt] = fs->make<TGraph>(mNGraphPoints, veta[ipt], vUnc_pt[ipt]);
0194 sprintf(name, "UNC_vs_Eta_CorPt%1.1f", mVPt[ipt]);
0195 mUncEta[ipt]->SetName(name);
0196 }
0197 for (unsigned ieta = 0; ieta < mVEta.size(); ieta++) {
0198 mVGraphPt[ieta] =
0199 fs->make<TGraphErrors>(mNGraphPoints, vpt[ieta], vjec_eta[ieta], vex_eta[ieta], vjecUnc_eta[ieta]);
0200 sprintf(name, "JEC_vs_RawPt_eta%1.1f", mVEta[ieta]);
0201 mVGraphPt[ieta]->SetName(name);
0202 mVGraphCorPt[ieta] =
0203 fs->make<TGraphErrors>(mNGraphPoints, vptcor[ieta], vjec_eta[ieta], vex_eta[ieta], vjecUnc_eta[ieta]);
0204 sprintf(name, "JEC_vs_CorPt_eta%1.1f", mVEta[ieta]);
0205 mVGraphCorPt[ieta]->SetName(name);
0206 mUncCorPt[ieta] = fs->make<TGraph>(mNGraphPoints, vptcor[ieta], vUnc_eta[ieta]);
0207 sprintf(name, "UNC_vs_CorPt_eta%1.1f", mVEta[ieta]);
0208 mUncCorPt[ieta]->SetName(name);
0209 }
0210 }
0211
0212
0213 DEFINE_FWK_MODULE(JetCorrectorDemo);