HiJetAlgos/plugins/UETableProducer.cc

0001 // system include files
0002 #include <iostream>
0003 #include <fstream>
0004 #include <vector>
0005 #include <string>
0006 #include <algorithm>
0007
0008 // user include files
0009 #include "FWCore/Framework/interface/Frameworkfwd.h"
0010 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0011
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/Framework/interface/EventSetup.h"
0014 #include "FWCore/Framework/interface/MakerMacros.h"
0015 #include "FWCore/Framework/interface/ESHandle.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017
0018 #include "FWCore/ServiceRegistry/interface/Service.h"
0019
0020 #include "CondFormats/DataRecord/interface/HeavyIonUERcd.h"
0021 #include "CondFormats/HIObjects/interface/UETable.h"
0022 #include "CondCore/DBOutputService/interface/PoolDBOutputService.h"
0023
0024 // For DB entry using JetCorrector to store the vector of float
0025 #include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
0026
0027 #include "FWCore/ParameterSet/interface/FileInPath.h"
0028
0029 using namespace std;
0030 //
0031 // class decleration
0032 //
0033
0034 namespace {
0035   class UETableProducer : public edm::one::EDAnalyzer<> {
0036   public:
0037     explicit UETableProducer(const edm::ParameterSet&);
0038     ~UETableProducer() override;
0039
0040     static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0041
0042   private:
0043     void beginJob() override {}
0044     void analyze(const edm::Event&, const edm::EventSetup&) override;
0045     void endJob() override;
0046
0047     // ----------member data ---------------------------
0048
0049     bool debug_;
0050     bool jetCorrectorFormat_;
0051
0052     string calibrationFile_;
0053     unsigned int runnum_;
0054
0055     unsigned int index = 0, np[5], ni0[2], ni1[2], ni2[2];
0056
0057     //    ue_interpolation_pf0[15][344],
0058     //    ue_interpolation_pf1[15][344],
0059     //    ue_interpolation_pf2[15][82];
0060   };
0061 }  // namespace
0062 //
0063 // constants, enums and typedefs
0064 //
0065
0066 //
0067 // static data member definitions
0068 //
0069
0070 //
0071 // constructors and destructor
0072 //
0073 UETableProducer::UETableProducer(const edm::ParameterSet& iConfig) : runnum_(0) {
0074   //now do what ever initialization is needed
0075   calibrationFile_ = iConfig.getParameter<std::string>("txtFile");
0076   //calibrationFile_ = "RecoHI/HiJetAlgos/data/ue_calibrations_pf_data.txt";
0077
0078   debug_ = iConfig.getUntrackedParameter<bool>("debug", false);
0079   jetCorrectorFormat_ = iConfig.getUntrackedParameter<bool>("jetCorrectorFormat", false);
0080
0081   np[0] = 3;   // Number of reduced PF ID (track, ECAL, HCAL)
0082   np[1] = 15;  // Number of pseudorapidity block
0083   np[2] = 5;   // Fourier series order
0084   np[3] = 2;   // Re or Im
0085   np[4] = 82;  // Number of feature parameter
0086
0087   ni0[0] = np[1];
0088   ni0[1] = 344;  // Number of track binning (kept = ECAL)
0089
0090   ni1[0] = np[1];
0091   ni1[1] = 344;  // Number of ECAL binning
0092
0093   ni2[0] = np[1];
0094   ni2[1] = 82;  // Number of HCAL binning
0095 }
0096
0097 UETableProducer::~UETableProducer() {
0098   // do anything here that needs to be done at desctruction time
0099   // (e.g. close files, deallocate resources etc.)
0100 }
0101
0102 //
0103 // member functions
0104 //
0105
0106 // ------------ method called to for each event  ------------
0107 void UETableProducer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0108   // nothing
0109 }
0110
0111 // ------------ method called once each job just after ending the event loop  ------------
0112 void UETableProducer::endJob() {
0113   std::string qpDataName = calibrationFile_;
0114   std::ifstream textTable_(qpDataName.c_str());
0115
0116   std::vector<float> ue_vec;
0117   std::unique_ptr<UETable> ue_predictor_pf;
0118
0119   if (!jetCorrectorFormat_) {
0120     ue_predictor_pf = std::make_unique<UETable>();
0121   }
0122   //  unsigned int Nnp_full = np[0] * np[1] * np[2] * np[3] * np[4];
0123   unsigned int Nnp = np[0] * np[1] * (1 + (np[2] - 1) * np[3]) * np[4];
0124   unsigned int Nni0 = ni0[0] * ni0[1];
0125   unsigned int Nni1 = ni1[0] * ni1[1];
0126   unsigned int Nni2 = ni2[0] * ni2[1];
0127
0128   if (!jetCorrectorFormat_) {
0129     ue_predictor_pf->np.resize(5);
0130     ue_predictor_pf->ni0.resize(2);
0131     ue_predictor_pf->ni1.resize(2);
0132     ue_predictor_pf->ni2.resize(2);
0133
0134     std::copy(np, np + 5, ue_predictor_pf->np.begin());
0135     std::copy(ni0, ni0 + 2, ue_predictor_pf->ni0.begin());
0136     std::copy(ni1, ni1 + 2, ue_predictor_pf->ni1.begin());
0137     std::copy(ni2, ni2 + 2, ue_predictor_pf->ni2.begin());
0138
0139     static const float edge_pseudorapidity[16] = {-5.191,
0140                                                   -2.650,
0141                                                   -2.043,
0142                                                   -1.740,
0143                                                   -1.479,
0144                                                   -1.131,
0145                                                   -0.783,
0146                                                   -0.522,
0147                                                   0.522,
0148                                                   0.783,
0149                                                   1.131,
0150                                                   1.479,
0151                                                   1.740,
0152                                                   2.043,
0153                                                   2.650,
0154                                                   5.191};
0155
0156     ue_predictor_pf->edgeEta.resize(16);
0157
0158     std::copy(edge_pseudorapidity, edge_pseudorapidity + 16, ue_predictor_pf->edgeEta.begin());
0159   }
0160
0161   std::string line;
0162
0163   while (std::getline(textTable_, line)) {
0164     if (line.empty() || line[0] == '#') {
0165       std::cout << " continue " << std::endl;
0166       continue;
0167     }
0168     std::istringstream linestream(line);
0169     float val;
0170     int bin0, bin1, bin2, bin3, bin4;
0171     if (index < Nnp) {
0172       linestream >> bin0 >> bin1 >> bin2 >> bin3 >> bin4 >> val;
0173       ue_vec.push_back(val);
0174     } else if (index < Nnp + Nni0) {
0175       linestream >> bin0 >> bin1 >> val;
0176       ue_vec.push_back(val);
0177     } else if (index < Nnp + Nni0 + Nni1) {
0178       linestream >> bin0 >> bin1 >> val;
0179       ue_vec.push_back(val);
0180     } else if (index < Nnp + Nni0 + Nni1 + Nni2) {
0181       linestream >> bin0 >> bin1 >> val;
0182       ue_vec.push_back(val);
0183     }
0184     ++index;
0185   }
0186
0187   edm::Service<cond::service::PoolDBOutputService> pool;
0188
0189   if (pool.isAvailable()) {
0190     if (jetCorrectorFormat_) {
0191       // A minimal dummy line that satisfies the JME # token >= 6 requirement, and has the correct key type
0192       JetCorrectorParameters::Definitions definition("1 0 0 0 Correction L1Offset");
0193       std::vector<JetCorrectorParameters::Record> record(
0194           1,
0195           JetCorrectorParameters::Record(
0196               ue_vec.size(), std::vector<float>(ue_vec.size(), 0), std::vector<float>(ue_vec.size(), 0), ue_vec));
0197       JetCorrectorParameters parameter(definition, record);
0198
0199       std::unique_ptr<JetCorrectorParametersCollection> jme_payload =
0200           std::make_unique<JetCorrectorParametersCollection>();
0201
0202       jme_payload->push_back(JetCorrectorParametersCollection::L1Offset, parameter);
0203
0204       if (pool->isNewTagRequest("JetCorrectionsRecord")) {
0205         pool->createOneIOV<JetCorrectorParametersCollection>(*jme_payload, pool->beginOfTime(), "JetCorrectionsRecord");
0206       } else {
0207         pool->appendOneIOV<JetCorrectorParametersCollection>(*jme_payload, pool->currentTime(), "JetCorrectionsRecord");
0208       }
0209     } else {
0210       ue_predictor_pf->values = ue_vec;
0211
0212       if (pool->isNewTagRequest("HeavyIonUERcd")) {
0213         pool->createOneIOV<UETable>(*ue_predictor_pf, pool->beginOfTime(), "HeavyIonUERcd");
0214       } else {
0215         pool->appendOneIOV<UETable>(*ue_predictor_pf, pool->currentTime(), "HeavyIonUERcd");
0216       }
0217     }
0218   }
0219 }
0220
0221 void UETableProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0222   edm::ParameterSetDescription desc;
0223   desc.add<std::string>("txtFile", "example");
0224   desc.addUntracked<bool>("debug", false);
0225   desc.addUntracked<bool>("jetCorrectorFormat", false);
0226   descriptions.add("produceUETable", desc);
0227 }
0228
0229 //define this as a plug-in
0230 DEFINE_FWK_MODULE(UETableProducer);