Hcal/bin/modGains.cc

0001 #include <cstdlib>
0002 #include <iostream>
0003 #include <fstream>
0004 #include <vector>
0005 #include <cstring>
0006 #include <string>
0007 #include <memory>
0008
0009 #include "CalibCalorimetry/HcalAlgos/interface/HcalDbASCIIIO.h"
0010 #include "CondFormats/HcalObjects/interface/HcalGains.h"
0011 #include "CondFormats/HcalObjects/interface/HcalRespCorrs.h"
0012 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0013 #include "FWCore/Framework/interface/Frameworkfwd.h"
0014 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0015 #include "FWCore/Framework/interface/Event.h"
0016 #include "FWCore/Framework/interface/MakerMacros.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 #include "FWCore/Utilities/interface/Exception.h"
0019 #include "Geometry/Records/interface/HcalRecNumberingRecord.h"
0020
0021 class modGains : public edm::one::EDAnalyzer<edm::one::WatchRuns> {
0022 public:
0023   explicit modGains(const edm::ParameterSet&);
0024   ~modGains() override;
0025
0026 private:
0027   void beginRun(edm::Run const& iEvent, edm::EventSetup const&) override {}
0028   void analyze(edm::Event const&, edm::EventSetup const&) override;
0029   void endRun(edm::Run const& iEvent, edm::EventSetup const&) override {}
0030   std::string s_operation;
0031   std::string fileIn, fileOut, fileCorr;
0032   double val;
0033   bool vectorop;
0034   edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> tok_htopo_;
0035 };
0036
0037 modGains::modGains(const edm::ParameterSet& iConfig) : vectorop(false) {
0038   s_operation = iConfig.getUntrackedParameter<std::string>("Operation");
0039   fileIn = iConfig.getUntrackedParameter<std::string>("FileIn");
0040   fileOut = iConfig.getUntrackedParameter<std::string>("FileOut");
0041   fileCorr = iConfig.getUntrackedParameter<std::string>("FileCorr");
0042   val = iConfig.getUntrackedParameter<double>("ScalarFactor");
0043   tok_htopo_ = esConsumes<HcalTopology, HcalRecNumberingRecord>();
0044
0045   if ((std::strcmp(s_operation.c_str(), "add") == 0) || (std::strcmp(s_operation.c_str(), "sub") == 0) ||
0046       (std::strcmp(s_operation.c_str(), "mult") == 0) ||
0047       (std::strcmp(s_operation.c_str(), "div") == 0)) {  // vector operation
0048     vectorop = true;
0049   } else if ((std::strcmp(s_operation.c_str(), "sadd") == 0) || (std::strcmp(s_operation.c_str(), "ssub") == 0) ||
0050              (std::strcmp(s_operation.c_str(), "smult") == 0) ||
0051              (std::strcmp(s_operation.c_str(), "sdiv") == 0)) {  // scalar operation
0052     std::cerr << "Scalar operation: using val=" << val << std::endl;
0053   } else {
0054     throw cms::Exception("Unknown", "modGains") << "Unknown operator. " << s_operation << " Stopping.\n";
0055   }
0056 }
0057
0058 modGains::~modGains() {}
0059
0060 void modGains::analyze(edm::Event const&, edm::EventSetup const& iSetup) {
0061   HcalTopology topo = iSetup.getData(tok_htopo_);
0062
0063   // get base conditions
0064   std::cerr << fileIn << std::endl;
0065   std::ifstream inStream(fileIn.c_str());
0066   HcalGains gainsIn(&topo);
0067   ;
0068   HcalDbASCIIIO::getObject(inStream, &gainsIn);
0069   inStream.close();
0070
0071   HcalRespCorrs corrsIn(&topo);
0072   ;
0073   if (vectorop) {
0074     std::ifstream inCorr(fileCorr.c_str());
0075     HcalDbASCIIIO::getObject(inCorr, &corrsIn);
0076     inCorr.close();
0077   }
0078
0079   HcalGains gainsOut(&topo);
0080   ;
0081   std::vector<DetId> channels = gainsIn.getAllChannels();
0082   std::cerr << "size = " << channels.size() << std::endl;
0083   for (unsigned i = 0; i < channels.size(); i++) {
0084     DetId id = channels[i];
0085
0086     if (vectorop) {  // vector operation
0087       if ((std::strcmp(s_operation.c_str(), "mult") == 0) || (std::strcmp(s_operation.c_str(), "div") == 0))
0088         val = 1.0;  // mult,div
0089       if ((std::strcmp(s_operation.c_str(), "add") == 0) || (std::strcmp(s_operation.c_str(), "sub") == 0))
0090         val = 0.0;  // add,sub
0091       if (corrsIn.exists(id)) {
0092         val = corrsIn.getValues(id)->getValue();
0093       }
0094       if (i % 100 == 0)
0095         std::cerr << "Vector operation, " << i << "th channel: using val=" << val << std::endl;
0096     }
0097
0098     std::unique_ptr<HcalGain> p_item;
0099     if ((std::strcmp(s_operation.c_str(), "add") == 0) || (std::strcmp(s_operation.c_str(), "sadd") == 0))
0100       p_item = std::make_unique<HcalGain>(id,
0101                                           gainsIn.getValues(id)->getValue(0) + val,
0102                                           gainsIn.getValues(id)->getValue(1) + val,
0103                                           gainsIn.getValues(id)->getValue(2) + val,
0104                                           gainsIn.getValues(id)->getValue(3) + val);
0105
0106     if ((std::strcmp(s_operation.c_str(), "sub") == 0) || (std::strcmp(s_operation.c_str(), "ssub") == 0))
0107       p_item = std::make_unique<HcalGain>(id,
0108                                           gainsIn.getValues(id)->getValue(0) - val,
0109                                           gainsIn.getValues(id)->getValue(1) - val,
0110                                           gainsIn.getValues(id)->getValue(2) - val,
0111                                           gainsIn.getValues(id)->getValue(3) - val);
0112
0113     if ((std::strcmp(s_operation.c_str(), "mult") == 0) || (std::strcmp(s_operation.c_str(), "smult") == 0))
0114       p_item = std::make_unique<HcalGain>(id,
0115                                           gainsIn.getValues(id)->getValue(0) * val,
0116                                           gainsIn.getValues(id)->getValue(1) * val,
0117                                           gainsIn.getValues(id)->getValue(2) * val,
0118                                           gainsIn.getValues(id)->getValue(3) * val);
0119
0120     if ((std::strcmp(s_operation.c_str(), "div") == 0) || (std::strcmp(s_operation.c_str(), "sdiv") == 0))
0121       p_item = std::make_unique<HcalGain>(id,
0122                                           gainsIn.getValues(id)->getValue(0) / val,
0123                                           gainsIn.getValues(id)->getValue(1) / val,
0124                                           gainsIn.getValues(id)->getValue(2) / val,
0125                                           gainsIn.getValues(id)->getValue(3) / val);
0126
0127     // for all
0128     if (p_item)
0129       gainsOut.addValues(*p_item);
0130     //    std::cerr << i << std::endl;
0131   }
0132   // write out
0133   std::ofstream outStream(fileOut.c_str());
0134   HcalDbASCIIIO::dumpObject(outStream, gainsOut);
0135   outStream.close();
0136 }
0137
0138 //define this as a plug-in
0139 DEFINE_FWK_MODULE(modGains);