EcalPedestalOffsets/src/EcalPedOffset.cc

0001 // TODO: fix header here?
0002 /**
0003  * \file EcalPedOffset.cc
0004  *
0005  * \author P. Govoni (pietro.govoni@cernNOSPAM.ch) - originally
0006  * \author S. Cooper (seth.cooper@cernNOSPAM.ch)
0007  * Last updated: @DATE@ @AUTHOR@
0008  *
0009  */
0010
0011 #include <fstream>
0012 #include <iostream>
0013 #include <memory>
0014
0015 #include "CalibCalorimetry/EcalPedestalOffsets/interface/EcalPedOffset.h"
0016
0017 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0018
0019 #include "OnlineDB/EcalCondDB/interface/EcalCondDBInterface.h"
0020 #include "OnlineDB/EcalCondDB/interface/MonPedestalOffsetsDat.h"
0021 #include "OnlineDB/EcalCondDB/interface/MonRunIOV.h"
0022 #include "OnlineDB/EcalCondDB/interface/RunIOV.h"
0023 #include "OnlineDB/EcalCondDB/interface/RunTag.h"
0024
0025 #include "DataFormats/EcalRawData/interface/EcalDCCHeaderBlock.h"
0026
0027 using namespace cms;
0028 using namespace edm;
0029
0030 //! ctor
0031 EcalPedOffset::EcalPedOffset(const ParameterSet &paramSet)
0032     : m_barrelDigiCollection(paramSet.getParameter<edm::InputTag>("EBdigiCollection")),
0033       m_endcapDigiCollection(paramSet.getParameter<edm::InputTag>("EEdigiCollection")),
0034       m_headerCollection(paramSet.getParameter<edm::InputTag>("headerCollection")),
0035       m_rawDataToken(consumes<EcalRawDataCollection>(m_headerCollection)),
0036       m_ebDigiToken(consumes<EBDigiCollection>(m_barrelDigiCollection)),
0037       m_eeDigiToken(consumes<EEDigiCollection>(m_endcapDigiCollection)),
0038       m_mappingToken(esConsumes()),
0039       m_xmlFile(paramSet.getParameter<std::string>("xmlFile")),
0040       m_DACmin(paramSet.getUntrackedParameter<int>("DACmin", 0)),
0041       m_DACmax(paramSet.getUntrackedParameter<int>("DACmax", 256)),
0042       m_RMSmax(paramSet.getUntrackedParameter<double>("RMSmax", 2)),
0043       m_bestPed(paramSet.getUntrackedParameter<int>("bestPed", 200)),
0044       m_dbHostName(paramSet.getUntrackedParameter<std::string>("dbHostName", "0")),
0045       m_dbName(paramSet.getUntrackedParameter<std::string>("dbName", "0")),
0046       m_dbUserName(paramSet.getUntrackedParameter<std::string>("dbUserName")),
0047       m_dbPassword(paramSet.getUntrackedParameter<std::string>("dbPassword")),
0048       m_dbHostPort(paramSet.getUntrackedParameter<int>("dbHostPort", 1521)),
0049       m_create_moniov(paramSet.getUntrackedParameter<bool>("createMonIOV", false)),
0050       m_location(paramSet.getUntrackedParameter<std::string>("location", "H4")),
0051       m_run(-1),
0052       m_plotting(paramSet.getParameter<std::string>("plotting")),
0053       m_maxSlopeAllowed_(paramSet.getUntrackedParameter<double>("maxSlopeAllowed", -29)),
0054       m_minSlopeAllowed_(paramSet.getUntrackedParameter<double>("minSlopeAllowed", -18)),
0055       m_maxChi2OverNDFAllowed_(paramSet.getUntrackedParameter<double>("maxChi2OverNDF", 5)) {
0056   edm::LogInfo("EcalPedOffset") << " reading "
0057                                 << " m_DACmin: " << m_DACmin << " m_DACmax: " << m_DACmax << " m_RMSmax: " << m_RMSmax
0058                                 << " m_bestPed: " << m_bestPed;
0059 }
0060
0061 // -----------------------------------------------------------------------------
0062
0063 //! dtor
0064 EcalPedOffset::~EcalPedOffset() {
0065   for (std::map<int, TPedValues *>::iterator mapIt = m_pedValues.begin(); mapIt != m_pedValues.end(); ++mapIt)
0066     delete mapIt->second;
0067   for (std::map<int, TPedResult *>::iterator mapIt = m_pedResult.begin(); mapIt != m_pedResult.end(); ++mapIt)
0068     delete mapIt->second;
0069 }
0070
0071 // -----------------------------------------------------------------------------
0072
0073 //! begin the run
0074 void EcalPedOffset::beginRun(Run const &, EventSetup const &eventSetup) {
0075   LogDebug("EcalPedOffset") << "entering beginRun...";
0076
0077   // get the electronics map
0078   const auto mappingHandle = eventSetup.getHandle(m_mappingToken);
0079   ecalElectronicsMap_ = mappingHandle.product();
0080 }
0081
0082 // -----------------------------------------------------------------------------
0083
0084 //! end the run
0085 void EcalPedOffset::endRun(Run const &, EventSetup const &eventSetup) {}
0086
0087 // -----------------------------------------------------------------------------
0088
0089 //! perform the analysis
0090 void EcalPedOffset::analyze(Event const &event, EventSetup const &eventSetup) {
0091   LogDebug("EcalPedOffset") << "entering analyze ...";
0092
0093   // get the headers
0094   // (one header for each supermodule)
0095   edm::Handle<EcalRawDataCollection> DCCHeaders;
0096   event.getByToken(m_rawDataToken, DCCHeaders);
0097
0098   std::map<int, int> DACvalues;
0099
0100   if (m_run == -1)
0101     m_run = event.id().run();
0102
0103   // loop over the headers
0104   for (EcalRawDataCollection::const_iterator headerItr = DCCHeaders->begin(); headerItr != DCCHeaders->end();
0105        ++headerItr) {
0106     EcalDCCHeaderBlock::EcalDCCEventSettings settings = headerItr->getEventSettings();
0107     int FEDid = 600 + headerItr->id();
0108     DACvalues[FEDid] = settings.ped_offset;
0109     LogDebug("EcalPedOffset") << "Found FED: " << FEDid << " in DCC header";
0110   }
0111
0112   bool barrelDigisFound = true;
0113   bool endcapDigisFound = true;
0114   // get the barrel digis
0115   // (one digi for each crystal)
0116   Handle<EBDigiCollection> barrelDigis;
0117   event.getByToken(m_ebDigiToken, barrelDigis);
0118   if (!barrelDigis.isValid()) {
0119     edm::LogError("EcalPedOffset") << "Error! can't get the product " << m_barrelDigiCollection
0120                                    << "; not reading barrel digis";
0121     barrelDigisFound = false;
0122   }
0123
0124   if (barrelDigis->empty()) {
0125     edm::LogInfo("EcalPedOffset") << "Size of EBDigiCollection is zero;"
0126                                   << " not reading barrel digis";
0127     barrelDigisFound = false;
0128   }
0129
0130   // get the endcap digis
0131   // (one digi for each crystal)
0132   Handle<EEDigiCollection> endcapDigis;
0133   event.getByToken(m_eeDigiToken, endcapDigis);
0134   if (!endcapDigis.isValid()) {
0135     edm::LogError("EcalPedOffset") << "Error! can't get the product " << m_endcapDigiCollection
0136                                    << "; not reading endcap digis";
0137     endcapDigisFound = false;
0138   }
0139
0140   if (endcapDigis->empty()) {
0141     edm::LogInfo("EcalPedOffset") << "Size of EEDigiCollection is zero;"
0142                                   << " not reading endcap digis";
0143     endcapDigisFound = false;
0144   }
0145
0146   if (barrelDigisFound)
0147     readDACs(barrelDigis, DACvalues);
0148   if (endcapDigisFound)
0149     readDACs(endcapDigis, DACvalues);
0150   if (!barrelDigisFound && !endcapDigisFound)
0151     edm::LogError("EcalPedOffset") << "No digis found in the event!";
0152 }
0153
0154 // -----------------------------------------------------------------------------
0155
0156 void EcalPedOffset::readDACs(const edm::Handle<EBDigiCollection> &pDigis, const std::map<int, int> &_DACvalues) {
0157   std::map<int, int> DACvalues = _DACvalues;
0158   // loop over the digis
0159   for (EBDigiCollection::const_iterator itdigi = pDigis->begin(); itdigi != pDigis->end(); ++itdigi) {
0160     int gainId = ((EBDataFrame)(*itdigi)).sample(0).gainId();
0161     EBDetId detId = EBDetId(itdigi->id());
0162     EcalElectronicsId elecId = ecalElectronicsMap_->getElectronicsId(detId);
0163     int FEDid = 600 + elecId.dccId();
0164     int crystalId = detId.ic();
0165
0166     // TODO: Behavior here
0167     if (DACvalues.find(FEDid) == DACvalues.end()) {
0168       edm::LogError("EcalPedOffset") << "Error! EB DCCid of digi does not "
0169                                         "match any DCCid found in DCC headers"
0170                                      << FEDid;
0171     }
0172
0173     if (!m_pedValues.count(FEDid)) {
0174       LogDebug("EcalPedOffset") << "Inserting new TPedValues object for FED:" << FEDid;
0175       m_pedValues[FEDid] = new TPedValues(m_RMSmax, m_bestPed);
0176     }
0177
0178     // loop over the samples
0179     for (int iSample = 0; iSample < EcalDataFrame::MAXSAMPLES; ++iSample) {
0180       m_pedValues[FEDid]->insert(
0181           gainId, crystalId, DACvalues[FEDid], ((EBDataFrame)(*itdigi)).sample(iSample).adc(), crystalId);
0182     }
0183
0184   }  // end loop over digis
0185 }
0186
0187 // -----------------------------------------------------------------------------
0188
0189 void EcalPedOffset::readDACs(const edm::Handle<EEDigiCollection> &pDigis, const std::map<int, int> &_DACvalues) {
0190   std::map<int, int> DACvalues = _DACvalues;
0191   // loop over the digis
0192   for (EEDigiCollection::const_iterator itdigi = pDigis->begin(); itdigi != pDigis->end(); ++itdigi) {
0193     int gainId = ((EEDataFrame)(*itdigi)).sample(0).gainId();
0194     // int gainId = itdigi->sample(0).gainId();
0195     EEDetId detId = EEDetId(itdigi->id());
0196     EcalElectronicsId elecId = ecalElectronicsMap_->getElectronicsId(detId);
0197     int FEDid = 600 + elecId.dccId();
0198     int crystalId = 25 * elecId.towerId() + 5 * elecId.stripId() + elecId.xtalId();
0199     int endcapCrystalId = 100 * elecId.towerId() + 5 * (elecId.stripId() - 1) + elecId.xtalId();
0200
0201     // TODO: Behavior here
0202     if (DACvalues.find(FEDid) == DACvalues.end()) {
0203       edm::LogError("EcalPedOffset") << "Error! EE DCCid of digi does not "
0204                                         "match any DCCid found in DCC headers: "
0205                                      << FEDid;
0206     }
0207
0208     if (!m_pedValues.count(FEDid)) {
0209       LogDebug("EcalPedOffset") << "Inserting new TPedValues object for FED:" << FEDid;
0210       m_pedValues[FEDid] = new TPedValues(m_RMSmax, m_bestPed);
0211     }
0212
0213     // loop over the samples
0214     for (int iSample = 0; iSample < EcalDataFrame::MAXSAMPLES; ++iSample) {
0215       m_pedValues[FEDid]->insert(
0216           gainId, crystalId, DACvalues[FEDid], ((EEDataFrame)(*itdigi)).sample(iSample).adc(), endcapCrystalId);
0217     }
0218
0219   }  // end loop over digis
0220 }
0221
0222 // -----------------------------------------------------------------------------
0223
0224 //! perform the minimization and write results
0225 void EcalPedOffset::endJob() {
0226   for (std::map<int, TPedValues *>::const_iterator smPeds = m_pedValues.begin(); smPeds != m_pedValues.end();
0227        ++smPeds) {
0228     m_pedResult[smPeds->first] = new TPedResult((smPeds->second)->terminate(m_DACmin, m_DACmax));
0229   }
0230   edm::LogInfo("EcalPedOffset") << " results map size " << m_pedResult.size();
0231   writeXMLFiles(m_xmlFile);
0232
0233   if (m_plotting != '0')
0234     makePlots();
0235   if (m_dbHostName != '0')
0236     writeDb();
0237 }
0238
0239 // -----------------------------------------------------------------------------
0240
0241 //! write the m_pedResult in the DB
0242 //! FIXME divide into sub-tasks
0243 void EcalPedOffset::writeDb() {
0244   LogDebug("EcalPedOffset") << " entering writeDb ...";
0245
0246   // connect to the database
0247   EcalCondDBInterface *DBconnection;
0248   try {
0249     LogInfo("EcalPedOffset") << "Opening DB connection with TNS_ADMIN ...";
0250     DBconnection = new EcalCondDBInterface(m_dbName, m_dbUserName, m_dbPassword);
0251   } catch (std::runtime_error &e) {
0252     LogError("EcalPedOffset") << e.what();
0253     if (!m_dbHostName.empty()) {
0254       try {
0255         LogInfo("EcalPedOffset") << "Opening DB connection without TNS_ADMIN ...";
0256         DBconnection = new EcalCondDBInterface(m_dbHostName, m_dbName, m_dbUserName, m_dbPassword, m_dbHostPort);
0257       } catch (std::runtime_error &e) {
0258         LogError("EcalPedOffset") << e.what();
0259         return;
0260       }
0261     } else
0262       return;
0263   }
0264
0265   // define the query for RunIOV to get the right place in the database
0266   RunTag runtag;
0267   LocationDef locdef;
0268   RunTypeDef rundef;
0269   locdef.setLocation(m_location);
0270
0271   runtag.setGeneralTag("PEDESTAL-OFFSET");
0272   rundef.setRunType("PEDESTAL-OFFSET");
0273   // rundef.setRunType ("TEST");
0274   // runtag.setGeneralTag ("TEST");
0275
0276   runtag.setLocationDef(locdef);
0277   runtag.setRunTypeDef(rundef);
0278
0279   run_t run = m_run;  // FIXME dal config file
0280   // RunIOV runiov = DBconnection->fetchRunIOV (&runtag, run);
0281   RunIOV runiov = DBconnection->fetchRunIOV(m_location, run);
0282
0283   // MonRunIOV
0284   MonVersionDef monverdef;
0285   monverdef.setMonitoringVersion("test01");
0286   MonRunTag montag;
0287   montag.setMonVersionDef(monverdef);
0288   montag.setGeneralTag("CMSSW");
0289
0290   subrun_t subrun = 1;  // hardcoded!
0291
0292   MonRunIOV moniov;
0293
0294   try {
0295     runtag = runiov.getRunTag();
0296     moniov = DBconnection->fetchMonRunIOV(&runtag, &montag, run, subrun);
0297   } catch (std::runtime_error &e) {
0298     if (m_create_moniov) {
0299       // if not already in the DB create a new MonRunIOV
0300       Tm startSubRun;
0301       startSubRun.setToCurrentGMTime();
0302       // setup the MonIOV
0303       moniov.setRunIOV(runiov);
0304       moniov.setSubRunNumber(subrun);
0305       moniov.setSubRunStart(startSubRun);
0306       moniov.setMonRunTag(montag);
0307       LogDebug("EcalPedOffset") << " creating a new MonRunIOV";
0308     } else {
0309       edm::LogError("EcalPedOffset") << " no MonRunIOV existing in the DB";
0310       edm::LogError("EcalPedOffset") << " the result will not be stored into the DB";
0311       if (DBconnection) {
0312         delete DBconnection;
0313       }
0314       return;
0315     }
0316   }
0317
0318   // create the table to be filled and the map to be inserted
0319   EcalLogicID ecid;
0320   std::map<EcalLogicID, MonPedestalOffsetsDat> DBdataset;
0321   MonPedestalOffsetsDat DBtable;
0322
0323   // fill the table
0324
0325   // loop over the super-modules
0326   for (std::map<int, TPedResult *>::const_iterator result = m_pedResult.begin(); result != m_pedResult.end();
0327        ++result) {
0328     // loop over the crystals
0329     for (int xtal = 0; xtal < 1700; ++xtal) {
0330       DBtable.setDACG1(result->second->m_DACvalue[2][xtal]);
0331       DBtable.setDACG6(result->second->m_DACvalue[1][xtal]);
0332       DBtable.setDACG12(result->second->m_DACvalue[0][xtal]);
0333       DBtable.setTaskStatus(true);  // FIXME to be set correctly
0334
0335       // fill the table
0336       if (DBconnection) {
0337         try {
0338           int fedid = result->first;
0339           int eid = m_pedValues[fedid]->getCrystalNumber(xtal);
0340           // If eid is zero, that crystal was not present in digis
0341           if (eid == 0)
0342             continue;
0343
0344           if (fedid >= 601 && fedid <= 609) {
0345             // Add the FEDid part in for DB
0346             eid = eid + 10000 * (fedid - 600);
0347             ecid = DBconnection->getEcalLogicID("EE_elec_crystal_number", eid);
0348           } else if (fedid >= 610 && fedid <= 627) {
0349             ecid = DBconnection->getEcalLogicID("EB_crystal_number", fedid - 610 + 19, eid);
0350           } else if (fedid >= 628 && fedid <= 645) {
0351             ecid = DBconnection->getEcalLogicID("EB_crystal_number", fedid - 628 + 1, eid);
0352           } else if (fedid >= 646 && fedid <= 654) {
0353             // Add the FEDid part in for DB
0354             eid = eid + 10000 * (fedid - 600);
0355             ecid = DBconnection->getEcalLogicID("EE_elec_crystal_number", eid);
0356           } else
0357             LogError("EcalPedOffset") << "FEDid is out of range 601-654";
0358
0359           DBdataset[ecid] = DBtable;
0360         } catch (std::runtime_error &e) {
0361           edm::LogError("EcalPedOffset") << e.what();
0362         }
0363       }
0364     }  // loop over the crystals
0365   }  // loop over the super-modules
0366
0367   // insert the map of tables in the database
0368   if (DBconnection) {
0369     try {
0370       LogDebug("EcalPedOffset") << "Inserting dataset ... " << std::flush;
0371       if (!DBdataset.empty())
0372         DBconnection->insertDataSet(&DBdataset, &moniov);
0373       LogDebug("EcalPedOffset") << "done.";
0374     } catch (std::runtime_error &e) {
0375       edm::LogError("EcalPedOffset") << e.what();
0376     }
0377   }
0378
0379   if (DBconnection) {
0380     delete DBconnection;
0381   }
0382 }
0383
0384 // -----------------------------------------------------------------------------
0385
0386 //! write the m_pedResults to XML files
0387 void EcalPedOffset::writeXMLFiles(std::string fileName) {
0388   // loop over the super-modules
0389   for (std::map<int, TPedResult *>::const_iterator smRes = m_pedResult.begin(); smRes != m_pedResult.end(); ++smRes) {
0390     std::string thisSMFileName = fileName;
0391     // open the output stream
0392     thisSMFileName += "_";
0393     thisSMFileName += intToString(smRes->first);
0394     thisSMFileName += ".xml";
0395     std::ofstream xml_outfile;
0396     xml_outfile.open(thisSMFileName.c_str());
0397
0398     // write the header file
0399     xml_outfile << "<offsets>" << std::endl;
0400     xml_outfile << "<PEDESTAL_OFFSET_RELEASE VERSION_ID = \"SM1_VER1\"> \n";
0401     xml_outfile << "  <RELEASE_ID>RELEASE_1</RELEASE_ID>\n";
0402     xml_outfile << "  <SUPERMODULE>";
0403     xml_outfile << smRes->first;
0404     xml_outfile << "</SUPERMODULE>\n";
0405     xml_outfile << "  <TIME_STAMP> 070705 </TIME_STAMP>" << std::endl;
0406
0407     // loop over the crystals
0408     for (int xtal = 0; xtal < 1700; ++xtal) {
0409       int crystalNumber = m_pedValues[smRes->first]->getCrystalNumber(xtal);
0410       if (crystalNumber == 0)
0411         continue;
0412       xml_outfile << "  <PEDESTAL_OFFSET>\n";
0413       xml_outfile << "    <HIGH>" << ((smRes->second)->m_DACvalue)[0][xtal] << "</HIGH>\n";
0414       xml_outfile << "    <MED>" << ((smRes->second)->m_DACvalue)[1][xtal] << "</MED>\n";
0415       xml_outfile << "    <LOW>" << ((smRes->second)->m_DACvalue)[2][xtal] << "</LOW>\n";
0416       xml_outfile << "    <CRYSTAL> " << crystalNumber << " </CRYSTAL>\n";
0417       xml_outfile << "  </PEDESTAL_OFFSET>" << std::endl;
0418     }
0419
0420     // close the open tags
0421     xml_outfile << " </PEDESTAL_OFFSET_RELEASE>" << std::endl;
0422     xml_outfile << "</offsets>" << std::endl;
0423     xml_outfile.close();
0424   }  // loop over the super-modules
0425 }
0426
0427 // -----------------------------------------------------------------------------
0428
0429 //! create the plots of the DAC pedestal trend
0430 void EcalPedOffset::makePlots() {
0431   LogDebug("EcalPedOffset") << " entering makePlots ...";
0432
0433   edm::LogInfo("EcalPedOffset") << " map size: " << m_pedValues.size();
0434
0435   // create the ROOT file
0436   m_plotting += ".root";
0437
0438   TFile *rootFile = new TFile(m_plotting.c_str(), "RECREATE");
0439
0440   // loop over the supermodules
0441   for (std::map<int, TPedValues *>::const_iterator smPeds = m_pedValues.begin(); smPeds != m_pedValues.end();
0442        ++smPeds) {
0443     // make a folder in the ROOT file
0444     char folderName[120];
0445     sprintf(folderName, "FED%02d", smPeds->first);
0446     rootFile->mkdir(folderName);
0447     smPeds->second->makePlots(rootFile, folderName, m_maxSlopeAllowed_, m_minSlopeAllowed_, m_maxChi2OverNDFAllowed_);
0448   }
0449
0450   rootFile->Close();
0451   delete rootFile;
0452
0453   LogDebug("EcalPedOffset") << " DONE";
0454 }
0455
0456 // -----------------------------------------------------------------------------
0457
0458 // convert an int to a string
0459 std::string EcalPedOffset::intToString(int num) {
0460   // outputs the number into the string stream and then flushes
0461   // the buffer (makes sure the output is put into the stream)
0462   std::ostringstream myStream;
0463   myStream << num << std::flush;
0464   return (myStream.str());  // returns the string form of the stringstream object
0465 }