Hcal/plugins/HcalGainsCheck.cc

0001 #include "CondTools/Hcal/interface/HcalGainsCheck.h"
0002
0003 HcalGainsCheck::HcalGainsCheck(edm::ParameterSet const& ps) {
0004   rootfile = ps.getUntrackedParameter<std::string>("rootfile", "null");
0005   outfile = ps.getUntrackedParameter<std::string>("outFile", "null");
0006   dumpupdate = ps.getUntrackedParameter<std::string>("dumpUpdateGainsTo", "null");
0007   dumprefs = ps.getUntrackedParameter<std::string>("dumpRefGainsTo", "null");
0008   emapflag = ps.getUntrackedParameter<bool>("checkEmap", false);
0009   validategainsflag = ps.getUntrackedParameter<bool>("validateGains", false);
0010   epsilon = ps.getUntrackedParameter<double>("deltaG", 1000000);
0011   m_tok1 = esConsumes<HcalGains, HcalGainsRcd>(edm::ESInputTag("", "update"));
0012   m_tok2 = esConsumes<HcalGains, HcalGainsRcd>(edm::ESInputTag("", "reference"));
0013   m_tokmap = esConsumes<HcalElectronicsMap, HcalElectronicsMapRcd>(edm::ESInputTag("", "reference"));
0014 }
0015
0016 void HcalGainsCheck::beginJob() {
0017   f = new TFile(rootfile.c_str(), "RECREATE");
0018
0019   //book histos:
0020   ocMapUp = new TH2F("ocMapUp", "occupancy_map_updated_gains", 83, -41.5, 41.5, 72, 0.5, 72.5);
0021   ocMapRef = new TH2F("ocMapRef", "occupancy_map_updated_gains", 83, -41.5, 41.5, 72, 0.5, 72.5);
0022   //  valMapUp;
0023   //  valMapRef;
0024
0025   diffUpRefCap0 = new TH1F("diffUpRefCap0", "difference_update_reference_Cap0", 100, -0.5, 0.5);
0026   ratioUpRefCap0 = new TH1F("ratioUpRefCap0", "ration_update_reference_Cap0", 100, 0.5, 1.5);
0027   gainsUpCap0 = new TH1F("gainsUpCap0", "gains_update_Cap0", 100, 0.0, 0.6);
0028   gainsRefCap0 = new TH1F("gainsRefCap0", "gains_reference_Cap0", 100, 0.0, 0.6);
0029
0030   diffUpRefCap1 = new TH1F("diffUpRefCap1", "difference_update_reference_Cap1", 100, -0.5, 0.5);
0031   ratioUpRefCap1 = new TH1F("ratioUpRefCap1", "ration_update_reference_Cap1", 100, 0.5, 1.5);
0032   gainsUpCap1 = new TH1F("gainsUpCap1", "gains_update_Cap1", 100, 0.0, 0.6);
0033   gainsRefCap1 = new TH1F("gainsRefCap1", "gains_reference_Cap1", 100, 0.0, 0.6);
0034
0035   diffUpRefCap2 = new TH1F("diffUpRefCap2", "difference_update_reference_Cap2", 100, -0.5, 0.5);
0036   ratioUpRefCap2 = new TH1F("ratioUpRefCap2", "ration_update_reference_Cap2", 100, 0.5, 1.5);
0037   gainsUpCap2 = new TH1F("gainsUpCap2", "gains_update_Cap2", 100, 0.0, 0.6);
0038   gainsRefCap2 = new TH1F("gainsRefCap2", "gains_reference_Cap2", 100, 0.0, 0.6);
0039
0040   diffUpRefCap3 = new TH1F("diffUpRefCap3", "difference_update_reference_Cap3", 100, -0.5, 0.5);
0041   ratioUpRefCap3 = new TH1F("ratioUpRefCap3", "ration_update_reference_Cap3", 100, 0.5, 1.5);
0042   gainsUpCap3 = new TH1F("gainsUpCap3", "gains_update_Cap3", 100, 0.0, 0.6);
0043   gainsRefCap3 = new TH1F("gainsRefCap3", "gains_reference_Cap3", 100, 0.0, 0.6);
0044
0045   //  gainsUpCap0vsEta = new TGraph("gainsUpCap0vsEta","gains_update_Cap0_vsEta",100,-41,0.6);
0046   //  gainsRefCap0vsEta = new TGraph("gainsRefCap0vsEta","gains_reference_Cap0_vsEta",100,0.0,0.6);
0047 }
0048
0049 void HcalGainsCheck::analyze(const edm::Event& ev, const edm::EventSetup& es) {
0050   using namespace edm::eventsetup;
0051   bool epsilonflag = false;
0052   bool notequalsflag = false;
0053   // get new gains
0054   const HcalGains* myNewGains = &es.getData(m_tok1);
0055
0056   // get reference gains
0057   const HcalGains* myRefGains = &es.getData(m_tok2);
0058
0059   // get e-map from reference
0060   const HcalElectronicsMap* myRefEMap = &es.getData(m_tokmap);
0061
0062   // dump gains:
0063   if (dumpupdate != "null") {
0064     std::ofstream outStream(dumpupdate.c_str());
0065     std::cout << "--- Dumping Gains - update ---" << std::endl;
0066     HcalDbASCIIIO::dumpObject(outStream, (*myNewGains));
0067   }
0068   if (dumprefs != "null") {
0069     std::ofstream outStream2(dumprefs.c_str());
0070     std::cout << "--- Dumping Gains - reference ---" << std::endl;
0071     HcalDbASCIIIO::dumpObject(outStream2, (*myRefGains));
0072   }
0073   // get the list of all channels
0074   std::vector<DetId> listNewChan = myNewGains->getAllChannels();
0075   std::vector<DetId> listRefChan = myRefGains->getAllChannels();
0076
0077   std::vector<DetId>::const_iterator cell;
0078
0079   //plots: occupancy map, value map, difference, ratio, gains:
0080   for (std::vector<DetId>::const_iterator it = listRefChan.begin(); it != listRefChan.end(); it++) {
0081     HcalGenericDetId myId(*it);
0082     //  ocMapRef->Fill(myId->);
0083
0084     float valCap0 = myRefGains->getValues(*it)->getValue(0);
0085     float valCap1 = myRefGains->getValues(*it)->getValue(1);
0086     float valCap2 = myRefGains->getValues(*it)->getValue(2);
0087     float valCap3 = myRefGains->getValues(*it)->getValue(3);
0088
0089     gainsRefCap0->Fill(valCap0);
0090     gainsRefCap1->Fill(valCap1);
0091     gainsRefCap2->Fill(valCap2);
0092     gainsRefCap3->Fill(valCap3);
0093
0094     cell = std::find(listNewChan.begin(), listNewChan.end(), (*it));
0095     if (cell != listNewChan.end())  //found
0096     {
0097       float valCap0up = myNewGains->getValues(*it)->getValue(0);
0098       float valCap1up = myNewGains->getValues(*it)->getValue(1);
0099       float valCap2up = myNewGains->getValues(*it)->getValue(2);
0100       float valCap3up = myNewGains->getValues(*it)->getValue(3);
0101
0102       diffUpRefCap0->Fill(valCap0up - valCap0);
0103       diffUpRefCap1->Fill(valCap1up - valCap1);
0104       diffUpRefCap2->Fill(valCap2up - valCap2);
0105       diffUpRefCap3->Fill(valCap3up - valCap3);
0106
0107       if (fabs(valCap0up - valCap0) > epsilon)
0108         epsilonflag = true;
0109       if (fabs(valCap1up - valCap1) > epsilon)
0110         epsilonflag = true;
0111       if (fabs(valCap2up - valCap2) > epsilon)
0112         epsilonflag = true;
0113       if (fabs(valCap3up - valCap3) > epsilon)
0114         epsilonflag = true;
0115
0116       if (valCap0up != valCap0)
0117         notequalsflag = true;
0118       if (valCap1up != valCap1)
0119         notequalsflag = true;
0120       if (valCap2up != valCap2)
0121         notequalsflag = true;
0122       if (valCap3up != valCap3)
0123         notequalsflag = true;
0124
0125       ratioUpRefCap0->Fill(valCap0up / valCap0);
0126       ratioUpRefCap1->Fill(valCap1up / valCap1);
0127       ratioUpRefCap2->Fill(valCap2up / valCap2);
0128       ratioUpRefCap3->Fill(valCap3up / valCap3);
0129     }
0130   }
0131   for (std::vector<DetId>::const_iterator it = listNewChan.begin(); it != listNewChan.end(); it++) {
0132     float valCap0 = myNewGains->getValues(*it)->getValue(0);
0133     float valCap1 = myNewGains->getValues(*it)->getValue(1);
0134     float valCap2 = myNewGains->getValues(*it)->getValue(2);
0135     float valCap3 = myNewGains->getValues(*it)->getValue(3);
0136
0137     gainsUpCap0->Fill(valCap0);
0138     gainsUpCap1->Fill(valCap1);
0139     gainsUpCap2->Fill(valCap2);
0140     gainsUpCap3->Fill(valCap3);
0141   }
0142
0143   if (epsilon != 1000000) {
0144     if (epsilonflag)
0145       throw cms::Exception("DataDoesNotMatch") << "Values differ by more than deltaG" << std::endl;
0146   } else {
0147     std::cout << "These gains do not differ by more than deltaG" << std::endl;
0148   }
0149
0150   if (validategainsflag) {
0151     if (notequalsflag)
0152       throw cms::Exception("DataDoesNotMatch") << "Values do not match" << std::endl;
0153   } else {
0154     std::cout << "These gains are identical" << std::endl;
0155   }
0156
0157   // go through list of valid channels from reference, look up if conditions exist for update
0158   // push back into new vector the corresponding updated conditions,
0159   // or if it doesn't exist, the reference
0160
0161   if (outfile != "null") {
0162     HcalGains* resultGains = new HcalGains(myRefGains->topo());
0163     for (std::vector<DetId>::const_iterator it = listRefChan.begin(); it != listRefChan.end(); it++) {
0164       DetId mydetid = *it;
0165       HcalGenericDetId myId(*it);
0166       cell = std::find(listNewChan.begin(), listNewChan.end(), mydetid);
0167       if (cell == listNewChan.end())  // not present in new list, take old conditions
0168       {
0169         const HcalGain* item = myRefGains->getValues(mydetid);
0170         std::cout << "o";
0171         resultGains->addValues(*item);
0172       } else  // present in new list, take new conditions
0173       {
0174         const HcalGain* item = myNewGains->getValues(mydetid);
0175         std::cout << "n";
0176         resultGains->addValues(*item);
0177       }
0178     }
0179     std::cout << std::endl;
0180
0181     std::vector<DetId> listResult = resultGains->getAllChannels();
0182     // get the e-map list of channels
0183     if (emapflag) {
0184       std::vector<HcalGenericDetId> listEMap = myRefEMap->allPrecisionId();
0185       // look up if emap channels are all present in pedestals, if not then cerr
0186       for (std::vector<HcalGenericDetId>::const_iterator it = listEMap.begin(); it != listEMap.end(); it++) {
0187         DetId mydetid = DetId(it->rawId());
0188         if (std::find(listResult.begin(), listResult.end(), mydetid) == listResult.end()) {
0189           std::cout << "Conditions not found for DetId = " << HcalGenericDetId(it->rawId()) << std::endl;
0190         }
0191       }
0192     }
0193
0194     // dump the resulting list of pedestals into a file
0195     //    std::ostringstream filename3;
0196     //    filename3 << "test_combined.txt";
0197     std::ofstream outStream3(outfile.c_str());
0198     std::cout << "--- Dumping Gains - the combined ones ---" << std::endl;
0199     HcalDbASCIIIO::dumpObject(outStream3, (*resultGains));
0200   }
0201   // const float* values = myped->getValues (channelID);
0202   //    if (values) std::cout << "pedestals for channel " << channelID << ": "
0203   //              << values [0] << '/' << values [1] << '/' << values [2] << '/' << values [3] << std::endl;
0204 }
0205
0206 // ------------ method called once each job just after ending the event loop  ------------
0207 void HcalGainsCheck::endJob() {
0208   if (rootfile != "null") {
0209     f->Write();
0210   }
0211   f->Close();
0212 }
0213
0214 DEFINE_FWK_MODULE(HcalGainsCheck);