Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:19:24 

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 "DataFormats/Candidate/interface/Particle.h"
0017 #include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h"
0018 #include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h"
0019 #include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
0020 #include "CondFormats/JetMETObjects/interface/FactorizedJetCorrector.h"
0021 //TFile Service
0022 #include "FWCore/ServiceRegistry/interface/Service.h"
0023 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0024 
0025 //
0026 // class declaration
0027 //
0028 
0029 class FactorizedJetCorrectorDemo : public edm::one::EDAnalyzer<> {
0030 public:
0031   explicit FactorizedJetCorrectorDemo(const edm::ParameterSet&);
0032   ~FactorizedJetCorrectorDemo() override;
0033   typedef reco::Particle::LorentzVector LorentzVector;
0034 
0035 private:
0036   void beginJob() override;
0037   void analyze(const edm::Event&, const edm::EventSetup&) override;
0038   void endJob() override;
0039 
0040   std::string mJetCorService, mPayloadName, mUncertaintyTag, mUncertaintyFile;
0041   edm::ESGetToken<JetCorrectorParametersCollection, JetCorrectionsRecord> mPayloadToken;
0042   std::vector<std::string> mLevels;
0043   bool mDebug, mUseCondDB;
0044   int mNHistoPoints, mNGraphPoints;
0045   double mEtaMin, mEtaMax, mPtMin, mPtMax;
0046   std::vector<double> mVEta, mVPt;
0047   double vjec_eta[100][1000], vjec_pt[100][1000], vpt[100][1000], vptcor[100][1000], veta[100][1000];
0048   double vjecUnc_eta[100][1000], vUnc_eta[100][1000], vjecUnc_pt[100][1000], vUnc_pt[100][1000], vex_eta[100][1000],
0049       vex_pt[100][1000];
0050   edm::Service<TFileService> fs;
0051   TH2F *mJECvsEta, *mJECvsPt;
0052   TGraphErrors *mVGraphEta[100], *mVGraphPt[100], *mVGraphCorPt[100];
0053   TGraph *mUncEta[100], *mUncCorPt[100];
0054   TRandom* mRandom;
0055 };
0056 //
0057 //----------- Class Implementation ------------------------------------------
0058 //
0059 //---------------------------------------------------------------------------
0060 FactorizedJetCorrectorDemo::FactorizedJetCorrectorDemo(const edm::ParameterSet& iConfig) {
0061   mLevels = iConfig.getParameter<std::vector<std::string> >("levels");
0062   mPayloadName = iConfig.getParameter<std::string>("PayloadName");
0063   mPayloadToken = esConsumes(edm::ESInputTag("", mPayloadName));
0064   mUncertaintyTag = iConfig.getParameter<std::string>("UncertaintyTag");
0065   mUncertaintyFile = iConfig.getParameter<std::string>("UncertaintyFile");
0066   mNHistoPoints = iConfig.getParameter<int>("NHistoPoints");
0067   mNGraphPoints = iConfig.getParameter<int>("NGraphPoints");
0068   mEtaMin = iConfig.getParameter<double>("EtaMin");
0069   mEtaMax = iConfig.getParameter<double>("EtaMax");
0070   mPtMin = iConfig.getParameter<double>("PtMin");
0071   mPtMax = iConfig.getParameter<double>("PtMax");
0072   mVEta = iConfig.getParameter<std::vector<double> >("VEta");
0073   mVPt = iConfig.getParameter<std::vector<double> >("VPt");
0074   mDebug = iConfig.getUntrackedParameter<bool>("Debug", false);
0075 }
0076 //---------------------------------------------------------------------------
0077 FactorizedJetCorrectorDemo::~FactorizedJetCorrectorDemo() {}
0078 //---------------------------------------------------------------------------
0079 void FactorizedJetCorrectorDemo::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0080   if (mDebug)
0081     std::cout << "Hello from FactorizedJetCorrectorDemo" << std::endl;
0082   // retreive parameters from the DB this still need a proper configurable
0083   // payloadName like: JetCorrectorParametersCollection_Spring10_AK5Calo.
0084   edm::ESHandle<JetCorrectorParametersCollection> parameters = iSetup.getHandle(mPayloadToken);
0085 
0086   std::vector<JetCorrectorParameters> params;
0087   for (std::vector<std::string>::const_iterator level = mLevels.begin(); level != mLevels.end(); ++level) {
0088     const JetCorrectorParameters& ip = (*parameters)[*level];  //ip.printScreen();
0089     if (mDebug)
0090       std::cout << "Adding level " << *level << std::endl;
0091     params.push_back(ip);
0092   }
0093 
0094   std::shared_ptr<FactorizedJetCorrector> corrector(new FactorizedJetCorrector(params));
0095 
0096   double jec, rawPt, corPt, eta;
0097   TLorentzVector P4;
0098   double dEta = (mEtaMax - mEtaMin) / mNGraphPoints;
0099   if (mDebug)
0100     std::cout << "Making JEC vs Eta and pT" << std::endl;
0101   for (int i = 0; i < mNHistoPoints; i++) {
0102     rawPt = mRandom->Uniform(mPtMin, mPtMax);
0103     eta = mRandom->Uniform(mEtaMin, mEtaMax);
0104     P4.SetPtEtaPhiE(rawPt, eta, 0, 0);
0105     corrector->setJetEta(eta);
0106     corrector->setJetPt(rawPt);
0107     corrector->setJetE(P4.E());
0108     jec = corrector->getCorrection();
0109     mJECvsEta->Fill(eta, jec);
0110     mJECvsPt->Fill(rawPt, jec);
0111   }
0112   if (mDebug)
0113     std::cout << "Making JEC vs pT for different etas" << std::endl;
0114   //--------- Pt Graphs ------------------
0115   for (unsigned ieta = 0; ieta < mVEta.size(); ieta++) {
0116     double rPt = pow((3500. / TMath::CosH(mVEta[ieta])) / mPtMin, 1. / mNGraphPoints);
0117     for (int i = 0; i < mNGraphPoints; i++) {
0118       rawPt = mPtMin * pow(rPt, i);
0119       eta = mVEta[ieta];
0120       vpt[ieta][i] = rawPt;
0121       P4.SetPtEtaPhiE(rawPt, eta, 0, 0);
0122       corrector->setJetEta(eta);
0123       corrector->setJetPt(rawPt);
0124       corrector->setJetE(P4.E());
0125       jec = corrector->getCorrection();
0126       vjec_eta[ieta][i] = jec;
0127       vptcor[ieta][i] = rawPt * jec;
0128       vex_eta[ieta][i] = 0.0;
0129       if (mDebug)
0130         std::cout << rawPt << " " << eta << " " << jec << " " << rawPt * jec << std::endl;
0131     }
0132   }
0133   if (mDebug)
0134     std::cout << "Making JEC vs eta for different pTs" << std::endl;
0135   //--------- Eta Graphs -------------------
0136   for (unsigned ipt = 0; ipt < mVPt.size(); ipt++) {
0137     for (int i = 0; i < mNGraphPoints; i++) {
0138       eta = mEtaMin + i * dEta;
0139       corPt = mVPt[ipt];
0140       veta[ipt][i] = eta;
0141       //---------- find the raw pt -----------
0142       double e = 1.0;
0143       int nLoop(0);
0144       rawPt = corPt;
0145       while (e > 0.0001 && nLoop < 10) {
0146         P4.SetPtEtaPhiE(rawPt, eta, 0, 0);
0147         LorentzVector rawP4(P4.Px(), P4.Py(), P4.Pz(), P4.E());
0148         corrector->setJetEta(eta);
0149         corrector->setJetPt(rawPt);
0150         corrector->setJetE(P4.E());
0151         jec = corrector->getCorrection();
0152         double tmp = rawPt * jec;
0153         e = fabs(tmp - corPt) / corPt;
0154         if (jec > 0)
0155           rawPt = corPt / jec;
0156         nLoop++;
0157       }
0158       //--------- calculate the jec for the rawPt --------
0159       P4.SetPtEtaPhiE(rawPt, eta, 0, 0);
0160       LorentzVector rawP4(P4.Px(), P4.Py(), P4.Pz(), P4.E());
0161       corrector->setJetEta(eta);
0162       corrector->setJetPt(rawPt);
0163       corrector->setJetE(P4.E());
0164       jec = corrector->getCorrection();
0165       vjec_pt[ipt][i] = jec;
0166       if (mDebug)
0167         std::cout << rawPt << " " << eta << " " << jec << " " << rawPt * jec << std::endl;
0168     }
0169   }
0170   if (mDebug)
0171     std::cout << "See ya!" << std::endl;
0172 }
0173 //---------------------------------------------------------------------------
0174 void FactorizedJetCorrectorDemo::beginJob() {
0175   if (mNGraphPoints > 1000)
0176     throw cms::Exception("FactorizedJetCorrectorDemo", "Too many graph points !!! Maximum is 1000 !!!");
0177   if (mVEta.size() > 100)
0178     throw cms::Exception("FactorizedJetCorrectorDemo", "Too many eta values !!! Maximum is 100 !!!");
0179   if (mVPt.size() > 100)
0180     throw cms::Exception("FactorizedJetCorrectorDemo", "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 FactorizedJetCorrectorDemo::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(FactorizedJetCorrectorDemo);
[ 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!