CMS 3D CMS Logo

/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/CondTools/Hcal/plugins/HcalGainsCheck.cc

Go to the documentation of this file.
00001 #include "CondTools/Hcal/interface/HcalGainsCheck.h"
00002 
00003 HcalGainsCheck::HcalGainsCheck(edm::ParameterSet const& ps)
00004 {
00005   rootfile = ps.getUntrackedParameter<std::string>("rootfile","null");
00006   outfile = ps.getUntrackedParameter<std::string>("outFile","null");
00007   dumpupdate = ps.getUntrackedParameter<std::string>("dumpUpdateGainsTo","null");
00008   dumprefs = ps.getUntrackedParameter<std::string>("dumpRefGainsTo","null");
00009   emapflag = ps.getUntrackedParameter<bool>("checkEmap",false);
00010   validategainsflag = ps.getUntrackedParameter<bool>("validateGains",false);
00011   epsilon = ps.getUntrackedParameter<double>("deltaG",1000000);
00012 }
00013 
00014 void HcalGainsCheck::beginJob()
00015 {
00016   f = new TFile(rootfile.c_str(),"RECREATE");
00017 
00018   //book histos:
00019   ocMapUp = new TH2F("ocMapUp","occupancy_map_updated_gains",83,-41.5,41.5,72,0.5,72.5);
00020   ocMapRef = new TH2F("ocMapRef","occupancy_map_updated_gains",83,-41.5,41.5,72,0.5,72.5);
00021 //  valMapUp;
00022 //  valMapRef;
00023 
00024   diffUpRefCap0 = new TH1F("diffUpRefCap0","difference_update_reference_Cap0",100,-0.5,0.5);
00025   ratioUpRefCap0 = new TH1F("ratioUpRefCap0", "ration_update_reference_Cap0",100,0.5,1.5);
00026   gainsUpCap0 = new TH1F("gainsUpCap0","gains_update_Cap0",100,0.0,0.6);
00027   gainsRefCap0 = new TH1F("gainsRefCap0","gains_reference_Cap0",100,0.0,0.6);
00028 
00029   diffUpRefCap1 = new TH1F("diffUpRefCap1","difference_update_reference_Cap1",100,-0.5,0.5);
00030   ratioUpRefCap1 = new TH1F("ratioUpRefCap1", "ration_update_reference_Cap1",100,0.5,1.5);
00031   gainsUpCap1 = new TH1F("gainsUpCap1","gains_update_Cap1",100,0.0,0.6);
00032   gainsRefCap1 = new TH1F("gainsRefCap1","gains_reference_Cap1",100,0.0,0.6);
00033 
00034   diffUpRefCap2 = new TH1F("diffUpRefCap2","difference_update_reference_Cap2",100,-0.5,0.5);
00035   ratioUpRefCap2 = new TH1F("ratioUpRefCap2", "ration_update_reference_Cap2",100,0.5,1.5);
00036   gainsUpCap2 = new TH1F("gainsUpCap2","gains_update_Cap2",100,0.0,0.6);
00037   gainsRefCap2 = new TH1F("gainsRefCap2","gains_reference_Cap2",100,0.0,0.6);
00038 
00039   diffUpRefCap3 = new TH1F("diffUpRefCap3","difference_update_reference_Cap3",100,-0.5,0.5);
00040   ratioUpRefCap3 = new TH1F("ratioUpRefCap3", "ration_update_reference_Cap3",100,0.5,1.5);
00041   gainsUpCap3 = new TH1F("gainsUpCap3","gains_update_Cap3",100,0.0,0.6);
00042   gainsRefCap3 = new TH1F("gainsRefCap3","gains_reference_Cap3",100,0.0,0.6);
00043 
00044   //  gainsUpCap0vsEta = new TGraph("gainsUpCap0vsEta","gains_update_Cap0_vsEta",100,-41,0.6);
00045   //  gainsRefCap0vsEta = new TGraph("gainsRefCap0vsEta","gains_reference_Cap0_vsEta",100,0.0,0.6);
00046 }
00047 
00048 
00049 void HcalGainsCheck::analyze(const edm::Event& ev, const edm::EventSetup& es)
00050 {
00051   using namespace edm::eventsetup;
00052   bool epsilonflag = false;
00053   bool notequalsflag = false;
00054   // get new gains
00055   edm::ESHandle<HcalGains> newGains;
00056   es.get<HcalGainsRcd>().get("update",newGains);
00057   const HcalGains* myNewGains = newGains.product();
00058 
00059   // get reference gains
00060   edm::ESHandle<HcalGains> refGains;
00061   es.get<HcalGainsRcd>().get("reference",refGains);
00062   const HcalGains* myRefGains = refGains.product();
00063 
00064   // get e-map from reference
00065   edm::ESHandle<HcalElectronicsMap> refEMap;
00066   es.get<HcalElectronicsMapRcd>().get("reference",refEMap);
00067   const HcalElectronicsMap* myRefEMap = refEMap.product();
00068 
00069 
00070   // dump gains:
00071    if(dumpupdate.compare("null")!=0){
00072    std::ofstream outStream(dumpupdate.c_str());
00073    std::cout << "--- Dumping Gains - update ---" << std::endl;
00074    HcalDbASCIIIO::dumpObject (outStream, (*myNewGains) );
00075    }
00076    if(dumprefs.compare("null")!=0){
00077     std::ofstream outStream2(dumprefs.c_str());
00078     std::cout << "--- Dumping Gains - reference ---" << std::endl;
00079     HcalDbASCIIIO::dumpObject (outStream2, (*myRefGains) );
00080    }
00081     // get the list of all channels
00082     std::vector<DetId> listNewChan = myNewGains->getAllChannels();
00083     std::vector<DetId> listRefChan = myRefGains->getAllChannels();
00084     
00085     std::vector<DetId>::const_iterator cell;
00086 
00087     //plots: occupancy map, value map, difference, ratio, gains:
00088     for (std::vector<DetId>::const_iterator it = listRefChan.begin(); it!=listRefChan.end(); it++)
00089       {
00090         HcalGenericDetId myId(*it);
00091         //      ocMapRef->Fill(myId->);
00092 
00093         float valCap0 = myRefGains->getValues(*it)->getValue(0);
00094         float valCap1 = myRefGains->getValues(*it)->getValue(1);
00095         float valCap2 = myRefGains->getValues(*it)->getValue(2);
00096         float valCap3 = myRefGains->getValues(*it)->getValue(3);
00097 
00098         gainsRefCap0->Fill(valCap0);
00099         gainsRefCap1->Fill(valCap1);
00100         gainsRefCap2->Fill(valCap2);
00101         gainsRefCap3->Fill(valCap3);
00102 
00103         cell = std::find(listNewChan.begin(), listNewChan.end(), (*it));
00104         if (cell != listNewChan.end() ) //found
00105           {
00106             float valCap0up = myNewGains->getValues(*it)->getValue(0);
00107             float valCap1up = myNewGains->getValues(*it)->getValue(1);
00108             float valCap2up = myNewGains->getValues(*it)->getValue(2);
00109             float valCap3up = myNewGains->getValues(*it)->getValue(3);
00110             
00111             diffUpRefCap0->Fill(valCap0up - valCap0);
00112             diffUpRefCap1->Fill(valCap1up - valCap1);
00113             diffUpRefCap2->Fill(valCap2up - valCap2);
00114             diffUpRefCap3->Fill(valCap3up - valCap3);
00115 
00116             if(fabs(valCap0up - valCap0) > epsilon) epsilonflag = true;
00117             if(fabs(valCap1up - valCap1) > epsilon) epsilonflag = true;
00118             if(fabs(valCap2up - valCap2) > epsilon) epsilonflag = true;
00119             if(fabs(valCap3up - valCap3) > epsilon) epsilonflag = true;
00120 
00121             if(valCap0up != valCap0) notequalsflag = true;
00122             if(valCap1up != valCap1) notequalsflag = true;
00123             if(valCap2up != valCap2) notequalsflag = true;
00124             if(valCap3up != valCap3) notequalsflag = true;
00125 
00126             ratioUpRefCap0->Fill(valCap0up / valCap0);
00127             ratioUpRefCap1->Fill(valCap1up / valCap1);
00128             ratioUpRefCap2->Fill(valCap2up / valCap2);
00129             ratioUpRefCap3->Fill(valCap3up / valCap3);
00130           }
00131       }
00132     for (std::vector<DetId>::const_iterator it = listNewChan.begin(); it!=listNewChan.end(); it++)
00133       {
00134         float valCap0 = myNewGains->getValues(*it)->getValue(0);
00135         float valCap1 = myNewGains->getValues(*it)->getValue(1);
00136         float valCap2 = myNewGains->getValues(*it)->getValue(2);
00137         float valCap3 = myNewGains->getValues(*it)->getValue(3);
00138 
00139         gainsUpCap0->Fill(valCap0);
00140         gainsUpCap1->Fill(valCap1);
00141         gainsUpCap2->Fill(valCap2);
00142         gainsUpCap3->Fill(valCap3);
00143       }
00144     
00145     if(epsilon != 1000000){
00146        if(epsilonflag) throw cms::Exception("DataDoesNotMatch") << "Values differ by more than deltaG" << std::endl;
00147     }else{
00148        std::cout << "These gains do not differ by more than deltaG" << std::endl;
00149     }
00150 
00151     if(validategainsflag){
00152        if(notequalsflag) throw cms::Exception("DataDoesNotMatch") << "Values do not match" << std::endl;
00153     }else{
00154        std::cout << "These gains are identical" << std::endl;
00155     }
00156 
00157     // go through list of valid channels from reference, look up if conditions exist for update
00158     // push back into new vector the corresponding updated conditions,
00159     // or if it doesn't exist, the reference
00160 
00161     if(outfile.compare("null")!=0){
00162       HcalGains *resultGains = new HcalGains(refGains->topo());
00163     for (std::vector<DetId>::const_iterator it = listRefChan.begin(); it != listRefChan.end(); it++)
00164       {
00165         DetId mydetid = *it;
00166         HcalGenericDetId myId(*it);
00167         cell = std::find(listNewChan.begin(), listNewChan.end(), mydetid);
00168         if (cell == listNewChan.end()) // not present in new list, take old conditions
00169           {
00170             const HcalGain* item = myRefGains->getValues(mydetid);
00171             std::cout << "o";
00172             resultGains->addValues(*item);
00173           }
00174         else // present in new list, take new conditions
00175           {
00176             const HcalGain* item = myNewGains->getValues(mydetid);
00177             std::cout << "n";
00178             resultGains->addValues(*item);
00179           }
00180       }
00181     std::cout << std::endl;
00182     
00183     std::vector<DetId> listResult = resultGains->getAllChannels();
00184     // get the e-map list of channels
00185     if(emapflag){
00186     std::vector<HcalGenericDetId> listEMap = myRefEMap->allPrecisionId();
00187     // look up if emap channels are all present in pedestals, if not then cerr
00188     for (std::vector<HcalGenericDetId>::const_iterator it = listEMap.begin(); it != listEMap.end(); it++)
00189       {
00190       DetId mydetid = DetId(it->rawId());
00191         if (std::find(listResult.begin(), listResult.end(), mydetid ) == listResult.end()  )
00192           {
00193             std::cout << "Conditions not found for DetId = " << HcalGenericDetId(it->rawId()) << std::endl;
00194           }
00195       }
00196     }
00197    
00198     // dump the resulting list of pedestals into a file
00199     //    std::ostringstream filename3;
00200     //    filename3 << "test_combined.txt";
00201     std::ofstream outStream3(outfile.c_str());
00202     std::cout << "--- Dumping Gains - the combined ones ---" << std::endl;
00203     HcalDbASCIIIO::dumpObject (outStream3, (*resultGains) );
00204 
00205     }
00206     // const float* values = myped->getValues (channelID);
00207     //    if (values) std::cout << "pedestals for channel " << channelID << ": "
00208     //                    << values [0] << '/' << values [1] << '/' << values [2] << '/' << values [3] << std::endl; 
00209 
00210 }
00211 
00212 
00213 // ------------ method called once each job just after ending the event loop  ------------
00214 void 
00215 HcalGainsCheck::endJob() 
00216 {
00217   if(rootfile.compare("null")!=0){
00218   f->Write();
00219 }
00220   f->Close();
00221 
00222 }
00223 
00224 DEFINE_FWK_MODULE(HcalGainsCheck);