Back to home page

Project CMSSW displayed by LXR
 
 

    

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 // user include files
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 //TFile Service
0021 #include "FWCore/ServiceRegistry/interface/Service.h"
0022 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0023 
0024 //
0025 // class declaration
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 //----------- Class Implementation ------------------------------------------
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   //--------- Pt Graphs ------------------
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);  // the jec is a function of the raw pt
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);  // the uncertainty is a function of the corrected pt
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   //--------- Eta Graphs -------------------
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       //---------- find the raw pt -----------
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       //--------- calculate the jec for the rawPt --------
0155       P4.SetPtEtaPhiE(rawPt, eta, 0, 0);
0156       LorentzVector rawP4(P4.Px(), P4.Py(), P4.Pz(), P4.E());
0157       jec = corrector->correction(rawP4);  // the jec is a function of the raw pt
0158       vjec_pt[ipt][i] = jec;
0159       unc = 0.0;
0160       if (!mUncertaintyTag.empty()) {
0161         jecUnc->setJetEta(eta);
0162         jecUnc->setJetPt(corPt);  // the uncertainty is a function of the corrected pt
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 //define this as a plug-in
0213 DEFINE_FWK_MODULE(JetCorrectorDemo);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!