EcalPedOffset.cc

Go to the documentation of this file.

// TODO: fix header here?
#include <fstream>
#include <iostream>
#include <memory>

#include "CalibCalorimetry/EcalPedestalOffsets/interface/EcalPedOffset.h"

#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "OnlineDB/EcalCondDB/interface/EcalCondDBInterface.h"
#include "OnlineDB/EcalCondDB/interface/MonPedestalOffsetsDat.h"
#include "OnlineDB/EcalCondDB/interface/MonRunIOV.h"
#include "OnlineDB/EcalCondDB/interface/RunIOV.h"
#include "OnlineDB/EcalCondDB/interface/RunTag.h"

#include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"

#include "DataFormats/EcalRawData/interface/EcalDCCHeaderBlock.h"
#include "DataFormats/EcalRawData/interface/EcalRawDataCollections.h"

#include "FWCore/Framework/interface/ESHandle.h"
#include "Geometry/EcalMapping/interface/EcalElectronicsMapping.h"
#include "Geometry/EcalMapping/interface/EcalMappingRcd.h"

using namespace cms;
using namespace edm;

EcalPedOffset::EcalPedOffset(const ParameterSet &paramSet)
    : m_barrelDigiCollection(paramSet.getParameter<edm::InputTag>("EBdigiCollection")),
      m_endcapDigiCollection(paramSet.getParameter<edm::InputTag>("EEdigiCollection")),
      m_headerCollection(paramSet.getParameter<edm::InputTag>("headerCollection")),
      m_xmlFile(paramSet.getParameter<std::string>("xmlFile")),
      m_DACmin(paramSet.getUntrackedParameter<int>("DACmin", 0)),
      m_DACmax(paramSet.getUntrackedParameter<int>("DACmax", 256)),
      m_RMSmax(paramSet.getUntrackedParameter<double>("RMSmax", 2)),
      m_bestPed(paramSet.getUntrackedParameter<int>("bestPed", 200)),
      m_dbHostName(paramSet.getUntrackedParameter<std::string>("dbHostName", "0")),
      m_dbName(paramSet.getUntrackedParameter<std::string>("dbName", "0")),
      m_dbUserName(paramSet.getUntrackedParameter<std::string>("dbUserName")),
      m_dbPassword(paramSet.getUntrackedParameter<std::string>("dbPassword")),
      m_dbHostPort(paramSet.getUntrackedParameter<int>("dbHostPort", 1521)),
      m_create_moniov(paramSet.getUntrackedParameter<bool>("createMonIOV", false)),
      m_location(paramSet.getUntrackedParameter<std::string>("location", "H4")),
      m_run(-1),
      m_plotting(paramSet.getParameter<std::string>("plotting")),
      m_maxSlopeAllowed_(paramSet.getUntrackedParameter<double>("maxSlopeAllowed", -29)),
      m_minSlopeAllowed_(paramSet.getUntrackedParameter<double>("minSlopeAllowed", -18)),
      m_maxChi2OverNDFAllowed_(paramSet.getUntrackedParameter<double>("maxChi2OverNDF", 5)) {
  edm::LogInfo("EcalPedOffset") << " reading "
                                << " m_DACmin: " << m_DACmin << " m_DACmax: " << m_DACmax << " m_RMSmax: " << m_RMSmax
                                << " m_bestPed: " << m_bestPed;
}

// -----------------------------------------------------------------------------

EcalPedOffset::~EcalPedOffset() {
  for (std::map<int, TPedValues *>::iterator mapIt = m_pedValues.begin(); mapIt != m_pedValues.end(); ++mapIt)
    delete mapIt->second;
  for (std::map<int, TPedResult *>::iterator mapIt = m_pedResult.begin(); mapIt != m_pedResult.end(); ++mapIt)
    delete mapIt->second;
}

// -----------------------------------------------------------------------------

void EcalPedOffset::beginRun(Run const &, EventSetup const &eventSetup) {
  LogDebug("EcalPedOffset") << "entering beginRun...";

  edm::ESHandle<EcalElectronicsMapping> handle;
  eventSetup.get<EcalMappingRcd>().get(handle);
  ecalElectronicsMap_ = handle.product();
}

// -----------------------------------------------------------------------------

void EcalPedOffset::analyze(Event const &event, EventSetup const &eventSetup) {
  LogDebug("EcalPedOffset") << "entering analyze ...";

  // get the headers
  // (one header for each supermodule)
  edm::Handle<EcalRawDataCollection> DCCHeaders;
  event.getByLabel(m_headerCollection, DCCHeaders);

  std::map<int, int> DACvalues;

  if (m_run == -1)
    m_run = event.id().run();

  // loop over the headers
  for (EcalRawDataCollection::const_iterator headerItr = DCCHeaders->begin(); headerItr != DCCHeaders->end();
       ++headerItr) {
    EcalDCCHeaderBlock::EcalDCCEventSettings settings = headerItr->getEventSettings();
    int FEDid = 600 + headerItr->id();
    DACvalues[FEDid] = settings.ped_offset;
    LogDebug("EcalPedOffset") << "Found FED: " << FEDid << " in DCC header";
  }

  bool barrelDigisFound = true;
  bool endcapDigisFound = true;
  // get the barrel digis
  // (one digi for each crystal)
  Handle<EBDigiCollection> barrelDigis;
  event.getByLabel(m_barrelDigiCollection, barrelDigis);
  if (!barrelDigis.isValid()) {
    edm::LogError("EcalPedOffset") << "Error! can't get the product " << m_barrelDigiCollection
                                   << "; not reading barrel digis";
    barrelDigisFound = false;
  }

  if (barrelDigis->empty()) {
    edm::LogInfo("EcalPedOffset") << "Size of EBDigiCollection is zero;"
                                  << " not reading barrel digis";
    barrelDigisFound = false;
  }

  // get the endcap digis
  // (one digi for each crystal)
  Handle<EEDigiCollection> endcapDigis;
  event.getByLabel(m_endcapDigiCollection, endcapDigis);
  if (!endcapDigis.isValid()) {
    edm::LogError("EcalPedOffset") << "Error! can't get the product " << m_endcapDigiCollection
                                   << "; not reading endcap digis";
    endcapDigisFound = false;
  }

  if (endcapDigis->empty()) {
    edm::LogInfo("EcalPedOffset") << "Size of EEDigiCollection is zero;"
                                  << " not reading endcap digis";
    endcapDigisFound = false;
  }

  if (barrelDigisFound)
    readDACs(barrelDigis, DACvalues);
  if (endcapDigisFound)
    readDACs(endcapDigis, DACvalues);
  if (!barrelDigisFound && !endcapDigisFound)
    edm::LogError("EcalPedOffset") << "No digis found in the event!";
}

// -----------------------------------------------------------------------------

void EcalPedOffset::readDACs(const edm::Handle<EBDigiCollection> &pDigis, const std::map<int, int> &_DACvalues) {
  std::map<int, int> DACvalues = _DACvalues;
  // loop over the digis
  for (EBDigiCollection::const_iterator itdigi = pDigis->begin(); itdigi != pDigis->end(); ++itdigi) {
    int gainId = ((EBDataFrame)(*itdigi)).sample(0).gainId();
    EBDetId detId = EBDetId(itdigi->id());
    EcalElectronicsId elecId = ecalElectronicsMap_->getElectronicsId(detId);
    int FEDid = 600 + elecId.dccId();
    int crystalId = detId.ic();

    // TODO: Behavior here
    if (DACvalues.find(FEDid) == DACvalues.end()) {
      edm::LogError("EcalPedOffset") << "Error! EB DCCid of digi does not "
                                        "match any DCCid found in DCC headers"
                                     << FEDid;
    }

    if (!m_pedValues.count(FEDid)) {
      LogDebug("EcalPedOffset") << "Inserting new TPedValues object for FED:" << FEDid;
      m_pedValues[FEDid] = new TPedValues(m_RMSmax, m_bestPed);
    }

    // loop over the samples
    for (int iSample = 0; iSample < EcalDataFrame::MAXSAMPLES; ++iSample) {
      m_pedValues[FEDid]->insert(
          gainId, crystalId, DACvalues[FEDid], ((EBDataFrame)(*itdigi)).sample(iSample).adc(), crystalId);
    }

  }  // end loop over digis
}

// -----------------------------------------------------------------------------

void EcalPedOffset::readDACs(const edm::Handle<EEDigiCollection> &pDigis, const std::map<int, int> &_DACvalues) {
  std::map<int, int> DACvalues = _DACvalues;
  // loop over the digis
  for (EEDigiCollection::const_iterator itdigi = pDigis->begin(); itdigi != pDigis->end(); ++itdigi) {
    int gainId = ((EEDataFrame)(*itdigi)).sample(0).gainId();
    // int gainId = itdigi->sample(0).gainId();
    EEDetId detId = EEDetId(itdigi->id());
    EcalElectronicsId elecId = ecalElectronicsMap_->getElectronicsId(detId);
    int FEDid = 600 + elecId.dccId();
    int crystalId = 25 * elecId.towerId() + 5 * elecId.stripId() + elecId.xtalId();
    int endcapCrystalId = 100 * elecId.towerId() + 5 * (elecId.stripId() - 1) + elecId.xtalId();

    // TODO: Behavior here
    if (DACvalues.find(FEDid) == DACvalues.end()) {
      edm::LogError("EcalPedOffset") << "Error! EE DCCid of digi does not "
                                        "match any DCCid found in DCC headers: "
                                     << FEDid;
    }

    if (!m_pedValues.count(FEDid)) {
      LogDebug("EcalPedOffset") << "Inserting new TPedValues object for FED:" << FEDid;
      m_pedValues[FEDid] = new TPedValues(m_RMSmax, m_bestPed);
    }

    // loop over the samples
    for (int iSample = 0; iSample < EcalDataFrame::MAXSAMPLES; ++iSample) {
      m_pedValues[FEDid]->insert(
          gainId, crystalId, DACvalues[FEDid], ((EEDataFrame)(*itdigi)).sample(iSample).adc(), endcapCrystalId);
    }

  }  // end loop over digis
}

// -----------------------------------------------------------------------------

void EcalPedOffset::endJob() {
  for (std::map<int, TPedValues *>::const_iterator smPeds = m_pedValues.begin(); smPeds != m_pedValues.end();
       ++smPeds) {
    m_pedResult[smPeds->first] = new TPedResult((smPeds->second)->terminate(m_DACmin, m_DACmax));
  }
  edm::LogInfo("EcalPedOffset") << " results map size " << m_pedResult.size();
  writeXMLFiles(m_xmlFile);

  if (m_plotting != '0')
    makePlots();
  if (m_dbHostName != '0')
    writeDb();
}

// -----------------------------------------------------------------------------

void EcalPedOffset::writeDb() {
  LogDebug("EcalPedOffset") << " entering writeDb ...";

  // connect to the database
  EcalCondDBInterface *DBconnection;
  try {
    LogInfo("EcalPedOffset") << "Opening DB connection with TNS_ADMIN ...";
    DBconnection = new EcalCondDBInterface(m_dbName, m_dbUserName, m_dbPassword);
  } catch (std::runtime_error &e) {
    LogError("EcalPedOffset") << e.what();
    if (!m_dbHostName.empty()) {
      try {
        LogInfo("EcalPedOffset") << "Opening DB connection without TNS_ADMIN ...";
        DBconnection = new EcalCondDBInterface(m_dbHostName, m_dbName, m_dbUserName, m_dbPassword, m_dbHostPort);
      } catch (std::runtime_error &e) {
        LogError("EcalPedOffset") << e.what();
        return;
      }
    } else
      return;
  }

  // define the query for RunIOV to get the right place in the database
  RunTag runtag;
  LocationDef locdef;
  RunTypeDef rundef;
  locdef.setLocation(m_location);

  runtag.setGeneralTag("PEDESTAL-OFFSET");
  rundef.setRunType("PEDESTAL-OFFSET");
  // rundef.setRunType ("TEST");
  // runtag.setGeneralTag ("TEST");

  runtag.setLocationDef(locdef);
  runtag.setRunTypeDef(rundef);

  run_t run = m_run;  // FIXME dal config file
  // RunIOV runiov = DBconnection->fetchRunIOV (&runtag, run);
  RunIOV runiov = DBconnection->fetchRunIOV(m_location, run);

  // MonRunIOV
  MonVersionDef monverdef;
  monverdef.setMonitoringVersion("test01");
  MonRunTag montag;
  montag.setMonVersionDef(monverdef);
  montag.setGeneralTag("CMSSW");

  subrun_t subrun = 1;  // hardcoded!

  MonRunIOV moniov;

  try {
    runtag = runiov.getRunTag();
    moniov = DBconnection->fetchMonRunIOV(&runtag, &montag, run, subrun);
  } catch (std::runtime_error &e) {
    if (m_create_moniov) {
      // if not already in the DB create a new MonRunIOV
      Tm startSubRun;
      startSubRun.setToCurrentGMTime();
      // setup the MonIOV
      moniov.setRunIOV(runiov);
      moniov.setSubRunNumber(subrun);
      moniov.setSubRunStart(startSubRun);
      moniov.setMonRunTag(montag);
      LogDebug("EcalPedOffset") << " creating a new MonRunIOV";
    } else {
      edm::LogError("EcalPedOffset") << " no MonRunIOV existing in the DB";
      edm::LogError("EcalPedOffset") << " the result will not be stored into the DB";
      if (DBconnection) {
        delete DBconnection;
      }
      return;
    }
  }

  // create the table to be filled and the map to be inserted
  EcalLogicID ecid;
  std::map<EcalLogicID, MonPedestalOffsetsDat> DBdataset;
  MonPedestalOffsetsDat DBtable;

  // fill the table

  // loop over the super-modules
  for (std::map<int, TPedResult *>::const_iterator result = m_pedResult.begin(); result != m_pedResult.end();
       ++result) {
    // loop over the crystals
    for (int xtal = 0; xtal < 1700; ++xtal) {
      DBtable.setDACG1(result->second->m_DACvalue[2][xtal]);
      DBtable.setDACG6(result->second->m_DACvalue[1][xtal]);
      DBtable.setDACG12(result->second->m_DACvalue[0][xtal]);
      DBtable.setTaskStatus(true);  // FIXME to be set correctly

      // fill the table
      if (DBconnection) {
        try {
          int fedid = result->first;
          int eid = m_pedValues[fedid]->getCrystalNumber(xtal);
          // If eid is zero, that crystal was not present in digis
          if (eid == 0)
            continue;

          if (fedid >= 601 && fedid <= 609) {
            // Add the FEDid part in for DB
            eid = eid + 10000 * (fedid - 600);
            ecid = DBconnection->getEcalLogicID("EE_elec_crystal_number", eid);
          } else if (fedid >= 610 && fedid <= 627) {
            ecid = DBconnection->getEcalLogicID("EB_crystal_number", fedid - 610 + 19, eid);
          } else if (fedid >= 628 && fedid <= 645) {
            ecid = DBconnection->getEcalLogicID("EB_crystal_number", fedid - 628 + 1, eid);
          } else if (fedid >= 646 && fedid <= 654) {
            // Add the FEDid part in for DB
            eid = eid + 10000 * (fedid - 600);
            ecid = DBconnection->getEcalLogicID("EE_elec_crystal_number", eid);
          } else
            LogError("EcalPedOffset") << "FEDid is out of range 601-654";

          DBdataset[ecid] = DBtable;
        } catch (std::runtime_error &e) {
          edm::LogError("EcalPedOffset") << e.what();
        }
      }
    }  // loop over the crystals
  }    // loop over the super-modules

  // insert the map of tables in the database
  if (DBconnection) {
    try {
      LogDebug("EcalPedOffset") << "Inserting dataset ... " << std::flush;
      if (!DBdataset.empty())
        DBconnection->insertDataSet(&DBdataset, &moniov);
      LogDebug("EcalPedOffset") << "done.";
    } catch (std::runtime_error &e) {
      edm::LogError("EcalPedOffset") << e.what();
    }
  }

  if (DBconnection) {
    delete DBconnection;
  }
}

// -----------------------------------------------------------------------------

void EcalPedOffset::writeXMLFiles(std::string fileName) {
  // loop over the super-modules
  for (std::map<int, TPedResult *>::const_iterator smRes = m_pedResult.begin(); smRes != m_pedResult.end(); ++smRes) {
    std::string thisSMFileName = fileName;
    // open the output stream
    thisSMFileName += "_";
    thisSMFileName += intToString(smRes->first);
    thisSMFileName += ".xml";
    std::ofstream xml_outfile;
    xml_outfile.open(thisSMFileName.c_str());

    // write the header file
    xml_outfile << "<offsets>" << std::endl;
    xml_outfile << "<PEDESTAL_OFFSET_RELEASE VERSION_ID = \"SM1_VER1\"> \n";
    xml_outfile << "  <RELEASE_ID>RELEASE_1</RELEASE_ID>\n";
    xml_outfile << "  <SUPERMODULE>";
    xml_outfile << smRes->first;
    xml_outfile << "</SUPERMODULE>\n";
    xml_outfile << "  <TIME_STAMP> 070705 </TIME_STAMP>" << std::endl;

    // loop over the crystals
    for (int xtal = 0; xtal < 1700; ++xtal) {
      int crystalNumber = m_pedValues[smRes->first]->getCrystalNumber(xtal);
      if (crystalNumber == 0)
        continue;
      xml_outfile << "  <PEDESTAL_OFFSET>\n";
      xml_outfile << "    <HIGH>" << ((smRes->second)->m_DACvalue)[0][xtal] << "</HIGH>\n";
      xml_outfile << "    <MED>" << ((smRes->second)->m_DACvalue)[1][xtal] << "</MED>\n";
      xml_outfile << "    <LOW>" << ((smRes->second)->m_DACvalue)[2][xtal] << "</LOW>\n";
      xml_outfile << "    <CRYSTAL> " << crystalNumber << " </CRYSTAL>\n";
      xml_outfile << "  </PEDESTAL_OFFSET>" << std::endl;
    }

    // close the open tags
    xml_outfile << " </PEDESTAL_OFFSET_RELEASE>" << std::endl;
    xml_outfile << "</offsets>" << std::endl;
    xml_outfile.close();
  }  // loop over the super-modules
}

// -----------------------------------------------------------------------------

void EcalPedOffset::makePlots() {
  LogDebug("EcalPedOffset") << " entering makePlots ...";

  edm::LogInfo("EcalPedOffset") << " map size: " << m_pedValues.size();

  // create the ROOT file
  m_plotting += ".root";

  TFile *rootFile = new TFile(m_plotting.c_str(), "RECREATE");

  // loop over the supermodules
  for (std::map<int, TPedValues *>::const_iterator smPeds = m_pedValues.begin(); smPeds != m_pedValues.end();
       ++smPeds) {
    // make a folder in the ROOT file
    char folderName[120];
    sprintf(folderName, "FED%02d", smPeds->first);
    rootFile->mkdir(folderName);
    smPeds->second->makePlots(rootFile, folderName, m_maxSlopeAllowed_, m_minSlopeAllowed_, m_maxChi2OverNDFAllowed_);
  }

  rootFile->Close();
  delete rootFile;

  LogDebug("EcalPedOffset") << " DONE";
}

// -----------------------------------------------------------------------------

// convert an int to a string
std::string EcalPedOffset::intToString(int num) {
  // outputs the number into the string stream and then flushes
  // the buffer (makes sure the output is put into the stream)
  std::ostringstream myStream;
  myStream << num << std::flush;
  return (myStream.str());  // returns the string form of the stringstream object
}