HcalCalibAlgos/test/HcalConstantsASCIIWriter.cc

0001 // system include files
0002 #include <fstream>
0003 #include <map>
0004 #include <iostream>
0005 #include <memory>
0006 #include <sstream>
0007 #include <string>
0008 #include <vector>
0009
0010 // user include files
0011 #include "FWCore/Framework/interface/Frameworkfwd.h"
0012 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/Framework/interface/EventSetup.h"
0015 #include "FWCore/Framework/interface/MakerMacros.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0017 #include "FWCore/ParameterSet/interface/FileInPath.h"
0018
0019 #include "CalibCalorimetry/HcalAlgos/interface/HcalDbASCIIIO.h"
0020 #include "CondFormats/HcalObjects/interface/HcalRespCorrs.h"
0021 #include "CondFormats/DataRecord/interface/HcalRespCorrsRcd.h"
0022
0023 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
0024 #include "DataFormats/Common/interface/Ref.h"
0025 #include "DataFormats/TrackReco/interface/Track.h"
0026 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
0027 #include "DataFormats/Candidate/interface/Candidate.h"
0028 #include "DataFormats/DetId/interface/DetId.h"
0029 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0030 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0031 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0032 #include "DataFormats/JetReco/interface/Jet.h"
0033 #include "DataFormats/JetReco/interface/CaloJet.h"
0034 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0035 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0036 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0037 #include "DataFormats/Provenance/interface/Provenance.h"
0038
0039 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0040 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0041 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0042 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0043 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0044 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0045
0046 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0047
0048 #include "HepMC/GenParticle.h"
0049 #include "HepMC/GenVertex.h"
0050
0051 //
0052 // class decleration
0053 //
0054 namespace cms {
0055   class HcalConstantsASCIIWriter : public edm::one::EDAnalyzer<> {
0056   public:
0057     explicit HcalConstantsASCIIWriter(const edm::ParameterSet&);
0058     ~HcalConstantsASCIIWriter();
0059
0060     virtual void analyze(const edm::Event&, const edm::EventSetup&);
0061     virtual void beginJob();
0062     virtual void endJob();
0063
0064   private:
0065     // ----------member data ---------------------------
0066     const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0067     const edm::ESGetToken<HcalRespCorrs, HcalRespCorrsRcd> tok_resp_;
0068
0069     std::ofstream* myout_hcal;
0070     std::string file_input;
0071     std::string file_output;
0072   };
0073 }  // namespace cms
0074
0075 //#define EDM_ML_DEBUG
0076
0077 //
0078 // constructors and destructor
0079 //
0080 namespace cms {
0081   HcalConstantsASCIIWriter::HcalConstantsASCIIWriter(const edm::ParameterSet& iConfig)
0082       : tok_geom_(esConsumes<CaloGeometry, CaloGeometryRecord>()),
0083         tok_resp_(esConsumes<HcalRespCorrs, HcalRespCorrsRcd>()) {
0084     // get name of output file with histogramms
0085     file_input = "Calibration/HcalCalibAlgos/data/" + iConfig.getParameter<std::string>("fileInput") + ".txt";
0086     file_output = "Calibration/HcalCalibAlgos/data/" + iConfig.getParameter<std::string>("fileOutput") + ".txt";
0087   }
0088
0089   HcalConstantsASCIIWriter::~HcalConstantsASCIIWriter() {
0090     // do anything here that needs to be done at desctruction time
0091     // (e.g. close files, deallocate resources etc.)
0092   }
0093
0094   void HcalConstantsASCIIWriter::beginJob() {
0095     edm::FileInPath f1(file_output);
0096     std::string fDataFile = f1.fullPath();
0097
0098     myout_hcal = new std::ofstream(fDataFile.c_str());
0099     if (!myout_hcal)
0100       edm::LogVerbatim("HcalCalib") << " Output file not open!!! ";
0101   }
0102
0103   void HcalConstantsASCIIWriter::endJob() { delete myout_hcal; }
0104
0105   //
0106   // member functions
0107   //
0108
0109   // ------------ method called to produce the data  ------------
0110   void HcalConstantsASCIIWriter::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0111     edm::LogVerbatim("HcalCalib") << " Start HcalConstantsASCIIWriter::analyze";
0112
0113     HcalRespCorrs* oldRespCorrs = new HcalRespCorrs(iSetup.getData(tok_resp_));
0114     //    std::vector<DetId> dd = oldRespCorrs->getAllChannels();
0115
0116     const CaloGeometry* geo = &iSetup.getData(tok_geom_);
0117     //   iSetup.get<HcalDbRecord>().get(conditions);
0118
0119     std::vector<DetId> did = geo->getValidDetIds();
0120
0121     std::map<HcalDetId, float> corrold;
0122     //map<HcalDetId,float> corrnew;
0123
0124     int mysubd, depth, ieta, iphi;
0125     float coradd, corerr;
0126
0127     std::vector<HcalDetId> theVector;
0128     for (std::vector<DetId>::iterator i = did.begin(); i != did.end(); i++) {
0129       if ((*i).det() == DetId::Hcal) {
0130         HcalDetId hid = HcalDetId(*i);
0131         theVector.push_back(hid);
0132         corrold[hid] = (oldRespCorrs->getValues(*i))->getValue();
0133         edm::LogVerbatim("HcalCalib") << " Old calibration " << hid.depth() << " " << hid.ieta() << " " << hid.iphi();
0134       }
0135     }
0136
0137     edm::LogVerbatim("HcalCalib") << " Get old calibration ";
0138     // Read new corrections from file
0139
0140     edm::FileInPath f1(file_input);
0141     std::string fDataFile = f1.fullPath();
0142
0143     std::ifstream in(fDataFile.c_str());
0144     std::string line;
0145
0146     double corrnew_p[5][5][45][75];
0147     double corrnew_m[5][5][45][75];
0148     edm::LogVerbatim("HcalCalib") << " Start to read txt file " << fDataFile.c_str() << std::endl;
0149     while (std::getline(in, line)) {
0150 #ifdef EDM_ML_DEBUG
0151       edm::LogVerbatim("HcalCalib") << " Line size " << line.size() << " " << line;
0152 #endif
0153
0154       if (!line.size() || line[0] == '#')
0155         continue;
0156       std::istringstream linestream(line);
0157
0158       linestream >> mysubd >> depth >> ieta >> iphi >> coradd >> corerr;
0159 #ifdef EDM_ML_DEBUG
0160       HcalDetId hid(HcalSubdetector(mysubd), ieta, iphi, depth);
0161       edm::LogVerbatim("HcalCalib") << " Check mysubd " << hid.subdet() << " depth " << hid.depth() << " ieta "
0162                                     << hid.ieta() << " iphi " << hid.iphi() << " " << hid.rawId();
0163 #endif
0164       int ietak = ieta;
0165       if (ieta < 0)
0166         ietak = -1 * ieta;
0167       if (ieta > 0)
0168         corrnew_p[mysubd][depth][ietak][iphi] = coradd;
0169       if (ieta < 0)
0170         corrnew_m[mysubd][depth][ietak][iphi] = coradd;
0171       edm::LogVerbatim("HcalCalib") << " Try to initialize mysubd " << mysubd << " depth " << depth << " ieta " << ieta
0172                                     << " " << ietak << " iphi " << iphi << " " << coradd;
0173     }
0174
0175     HcalRespCorrs* mycorrections = new HcalRespCorrs(oldRespCorrs->topo());
0176
0177     for (std::vector<HcalDetId>::iterator it = theVector.begin(); it != theVector.end(); it++) {
0178       float cc1 = (*corrold.find(*it)).second;
0179       //    float cc2 = (*corrnew.find(*it)).second;
0180       float cc2 = 0.;
0181       int ietak = (*it).ieta();
0182
0183       if ((*it).ieta() < 0)
0184         ietak = -1 * (*it).ieta();
0185
0186       if ((*it).ieta() > 0)
0187         cc2 = corrnew_p[(*it).subdet()][(*it).depth()][ietak][(*it).iphi()];
0188       if ((*it).ieta() < 0)
0189         cc2 = corrnew_m[(*it).subdet()][(*it).depth()][ietak][(*it).iphi()];
0190
0191       float cc = cc1 * cc2;
0192       edm::LogVerbatim("HcalCalib") << " Multiply " << (*it).subdet() << " " << (*it).depth() << " " << (*it).ieta()
0193                                     << " " << ietak << " " << (*it).iphi() << " " << (*it).rawId() << " " << cc1 << " "
0194                                     << cc2;
0195
0196       // now make the basic object for one cell with HcalDetId myDetId containing the value myValue
0197       HcalRespCorr item((*it).rawId(), cc);
0198       mycorrections->addValues(item);
0199     }
0200
0201     HcalRespCorrs mycc = *mycorrections;
0202     HcalDbASCIIIO::dumpObject(*myout_hcal, mycc);
0203   }
0204 }  // namespace cms
0205 //define this as a plug-in
0206
0207 #include "FWCore/Framework/interface/MakerMacros.h"
0208
0209 using cms::HcalConstantsASCIIWriter;
0210 DEFINE_FWK_MODULE(HcalConstantsASCIIWriter);