EcalRecHitProducer.cc

Go to the documentation of this file.

#include "RecoLocalCalo/EcalRecProducers/plugins/EcalRecHitProducer.h"
#include "RecoLocalCalo/EcalRecAlgos/interface/EcalCleaningAlgo.h"
#include "DataFormats/Common/interface/Handle.h"
#include "FWCore/Framework/interface/ESHandle.h"

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

#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
#include "DataFormats/EcalDetId/interface/EBDetId.h"
#include "DataFormats/EcalDetId/interface/EEDetId.h"
#include "DataFormats/EcalDetId/interface/EcalElectronicsId.h"
#include "DataFormats/EcalDetId/interface/EcalTrigTowerDetId.h"
#include "DataFormats/EcalDetId/interface/EcalScDetId.h"

#include "CondFormats/DataRecord/interface/EcalChannelStatusRcd.h"
#include "CondFormats/EcalObjects/interface/EcalChannelStatus.h"

#include "RecoLocalCalo/EcalRecProducers/interface/EcalRecHitWorkerFactory.h"

#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
#include "FWCore/ParameterSet/interface/ParameterSetDescription.h"

EcalRecHitProducer::EcalRecHitProducer(const edm::ParameterSet& ps) {
  ebRechitCollection_ = ps.getParameter<std::string>("EBrechitCollection");
  eeRechitCollection_ = ps.getParameter<std::string>("EErechitCollection");

  recoverEBIsolatedChannels_ = ps.getParameter<bool>("recoverEBIsolatedChannels");
  recoverEEIsolatedChannels_ = ps.getParameter<bool>("recoverEEIsolatedChannels");
  recoverEBVFE_ = ps.getParameter<bool>("recoverEBVFE");
  recoverEEVFE_ = ps.getParameter<bool>("recoverEEVFE");
  recoverEBFE_ = ps.getParameter<bool>("recoverEBFE");
  recoverEEFE_ = ps.getParameter<bool>("recoverEEFE");
  killDeadChannels_ = ps.getParameter<bool>("killDeadChannels");

  produces<EBRecHitCollection>(ebRechitCollection_);
  produces<EERecHitCollection>(eeRechitCollection_);

  ebUncalibRecHitToken_ =
      consumes<EBUncalibratedRecHitCollection>(ps.getParameter<edm::InputTag>("EBuncalibRecHitCollection"));

  eeUncalibRecHitToken_ =
      consumes<EEUncalibratedRecHitCollection>(ps.getParameter<edm::InputTag>("EEuncalibRecHitCollection"));

  ebDetIdToBeRecoveredToken_ = consumes<std::set<EBDetId>>(ps.getParameter<edm::InputTag>("ebDetIdToBeRecovered"));

  eeDetIdToBeRecoveredToken_ = consumes<std::set<EEDetId>>(ps.getParameter<edm::InputTag>("eeDetIdToBeRecovered"));

  ebFEToBeRecoveredToken_ = consumes<std::set<EcalTrigTowerDetId>>(ps.getParameter<edm::InputTag>("ebFEToBeRecovered"));

  eeFEToBeRecoveredToken_ = consumes<std::set<EcalScDetId>>(ps.getParameter<edm::InputTag>("eeFEToBeRecovered"));

  std::string componentType = ps.getParameter<std::string>("algo");
  edm::ConsumesCollector c{consumesCollector()};
  worker_ = std::unique_ptr<EcalRecHitWorkerBaseClass>{EcalRecHitWorkerFactory::get()->create(componentType, ps, c)};

  // to recover problematic channels
  componentType = ps.getParameter<std::string>("algoRecover");
  workerRecover_ = std::unique_ptr<EcalRecHitWorkerBaseClass>{EcalRecHitWorkerFactory::get()->create(componentType, ps, c)};

  edm::ParameterSet cleaningPs = ps.getParameter<edm::ParameterSet>("cleaningConfig");
  cleaningAlgo_ = std::make_unique<EcalCleaningAlgo>(cleaningPs);
}

EcalRecHitProducer::~EcalRecHitProducer() = default;

void EcalRecHitProducer::produce(edm::Event& evt, const edm::EventSetup& es) {
  using namespace edm;

  Handle<EBUncalibratedRecHitCollection> pEBUncalibRecHits;
  Handle<EEUncalibratedRecHitCollection> pEEUncalibRecHits;

  const EBUncalibratedRecHitCollection* ebUncalibRecHits = nullptr;
  const EEUncalibratedRecHitCollection* eeUncalibRecHits = nullptr;

  // get the barrel uncalib rechit collection

  evt.getByToken(ebUncalibRecHitToken_, pEBUncalibRecHits);
  ebUncalibRecHits = pEBUncalibRecHits.product();
  LogDebug("EcalRecHitDebug") << "total # EB uncalibrated rechits: " << ebUncalibRecHits->size();

  evt.getByToken(eeUncalibRecHitToken_, pEEUncalibRecHits);
  eeUncalibRecHits = pEEUncalibRecHits.product();  // get a ptr to the product
  LogDebug("EcalRecHitDebug") << "total # EE uncalibrated rechits: " << eeUncalibRecHits->size();

  // collection of rechits to put in the event
  auto ebRecHits = std::make_unique<EBRecHitCollection>();
  auto eeRecHits = std::make_unique<EERecHitCollection>();

  worker_->set(es);

  if (recoverEBIsolatedChannels_ || recoverEEIsolatedChannels_ || recoverEBFE_ || recoverEEFE_ || recoverEBVFE_ ||
      recoverEEVFE_ || killDeadChannels_) {
    workerRecover_->set(es);
  }

  if (ebUncalibRecHits) {
    // loop over uncalibrated rechits to make calibrated ones
    for (EBUncalibratedRecHitCollection::const_iterator it = ebUncalibRecHits->begin(); it != ebUncalibRecHits->end();
         ++it) {
      worker_->run(evt, *it, *ebRecHits);
    }
  }

  if (eeUncalibRecHits) {
    // loop over uncalibrated rechits to make calibrated ones
    for (EEUncalibratedRecHitCollection::const_iterator it = eeUncalibRecHits->begin(); it != eeUncalibRecHits->end();
         ++it) {
      worker_->run(evt, *it, *eeRecHits);
    }
  }

  // sort collections before attempting recovery, to avoid insertion of double recHits
  ebRecHits->sort();
  eeRecHits->sort();

  if (recoverEBIsolatedChannels_ || recoverEBFE_ || killDeadChannels_) {
    edm::Handle<std::set<EBDetId>> pEBDetId;
    const std::set<EBDetId>* detIds = nullptr;
    evt.getByToken(ebDetIdToBeRecoveredToken_, pEBDetId);
    detIds = pEBDetId.product();

    if (detIds) {
      edm::ESHandle<EcalChannelStatus> chStatus;
      es.get<EcalChannelStatusRcd>().get(chStatus);
      for (std::set<EBDetId>::const_iterator it = detIds->begin(); it != detIds->end(); ++it) {
        // get channel status map to treat dead VFE separately
        EcalChannelStatusMap::const_iterator chit = chStatus->find(*it);
        EcalChannelStatusCode chStatusCode;
        if (chit != chStatus->end()) {
          chStatusCode = *chit;
        } else {
          edm::LogError("EcalRecHitProducerError") << "No channel status found for xtal " << (*it).rawId()
                                                   << "! something wrong with EcalChannelStatus in your DB? ";
        }
        EcalUncalibratedRecHit urh;
        if (chStatusCode.getStatusCode() == EcalChannelStatusCode::kDeadVFE) {  // dead VFE (from DB info)
          // uses the EcalUncalibratedRecHit to pass the DetId info
          urh = EcalUncalibratedRecHit(*it, 0, 0, 0, 0, EcalRecHitWorkerBaseClass::EB_VFE);
          if (recoverEBVFE_ || killDeadChannels_)
            workerRecover_->run(evt, urh, *ebRecHits);
        } else {
          // uses the EcalUncalibratedRecHit to pass the DetId info
          urh = EcalUncalibratedRecHit(*it, 0, 0, 0, 0, EcalRecHitWorkerBaseClass::EB_single);
          if (recoverEBIsolatedChannels_ || killDeadChannels_)
            workerRecover_->run(evt, urh, *ebRecHits);
        }
      }
    }
  }

  if (recoverEEIsolatedChannels_ || recoverEEVFE_ || killDeadChannels_) {
    edm::Handle<std::set<EEDetId>> pEEDetId;
    const std::set<EEDetId>* detIds = nullptr;

    evt.getByToken(eeDetIdToBeRecoveredToken_, pEEDetId);
    detIds = pEEDetId.product();

    if (detIds) {
      edm::ESHandle<EcalChannelStatus> chStatus;
      es.get<EcalChannelStatusRcd>().get(chStatus);
      for (std::set<EEDetId>::const_iterator it = detIds->begin(); it != detIds->end(); ++it) {
        // get channel status map to treat dead VFE separately
        EcalChannelStatusMap::const_iterator chit = chStatus->find(*it);
        EcalChannelStatusCode chStatusCode;
        if (chit != chStatus->end()) {
          chStatusCode = *chit;
        } else {
          edm::LogError("EcalRecHitProducerError") << "No channel status found for xtal " << (*it).rawId()
                                                   << "! something wrong with EcalChannelStatus in your DB? ";
        }
        EcalUncalibratedRecHit urh;
        if (chStatusCode.getStatusCode() == EcalChannelStatusCode::kDeadVFE) {  // dead VFE (from DB info)
          // uses the EcalUncalibratedRecHit to pass the DetId info
          urh = EcalUncalibratedRecHit(*it, 0, 0, 0, 0, EcalRecHitWorkerBaseClass::EE_VFE);
          if (recoverEEVFE_ || killDeadChannels_)
            workerRecover_->run(evt, urh, *eeRecHits);
        } else {
          // uses the EcalUncalibratedRecHit to pass the DetId info
          urh = EcalUncalibratedRecHit(*it, 0, 0, 0, 0, EcalRecHitWorkerBaseClass::EE_single);
          if (recoverEEIsolatedChannels_ || killDeadChannels_)
            workerRecover_->run(evt, urh, *eeRecHits);
        }
      }
    }
  }

  if (recoverEBFE_ || killDeadChannels_) {
    edm::Handle<std::set<EcalTrigTowerDetId>> pEBFEId;
    const std::set<EcalTrigTowerDetId>* ttIds = nullptr;

    evt.getByToken(ebFEToBeRecoveredToken_, pEBFEId);
    ttIds = pEBFEId.product();

    if (ttIds) {
      for (std::set<EcalTrigTowerDetId>::const_iterator it = ttIds->begin(); it != ttIds->end(); ++it) {
        // uses the EcalUncalibratedRecHit to pass the DetId info
        int ieta = (((*it).ietaAbs() - 1) * 5 + 1) * (*it).zside();  // from EcalTrigTowerConstituentsMap
        int iphi = (((*it).iphi() - 1) * 5 + 11) % 360;              // from EcalTrigTowerConstituentsMap
        if (iphi <= 0)
          iphi += 360;  // from EcalTrigTowerConstituentsMap
        EcalUncalibratedRecHit urh(
            EBDetId(ieta, iphi, EBDetId::ETAPHIMODE), 0, 0, 0, 0, EcalRecHitWorkerBaseClass::EB_FE);
        workerRecover_->run(evt, urh, *ebRecHits);
      }
    }
  }

  if (recoverEEFE_ || killDeadChannels_) {
    edm::Handle<std::set<EcalScDetId>> pEEFEId;
    const std::set<EcalScDetId>* scIds = nullptr;

    evt.getByToken(eeFEToBeRecoveredToken_, pEEFEId);
    scIds = pEEFEId.product();

    if (scIds) {
      for (std::set<EcalScDetId>::const_iterator it = scIds->begin(); it != scIds->end(); ++it) {
        // uses the EcalUncalibratedRecHit to pass the DetId info
        if (EEDetId::validDetId(((*it).ix() - 1) * 5 + 1, ((*it).iy() - 1) * 5 + 1, (*it).zside())) {
          EcalUncalibratedRecHit urh(EEDetId(((*it).ix() - 1) * 5 + 1, ((*it).iy() - 1) * 5 + 1, (*it).zside()),
                                     0,
                                     0,
                                     0,
                                     0,
                                     EcalRecHitWorkerBaseClass::EE_FE);
          workerRecover_->run(evt, urh, *eeRecHits);
        }
      }
    }
  }

  // without re-sorting, find (used below in cleaning) will lead
  // to undefined results
  ebRecHits->sort();
  eeRecHits->sort();

  // apply spike cleaning
  if (cleaningAlgo_) {
    cleaningAlgo_->setFlags(*ebRecHits);
    cleaningAlgo_->setFlags(*eeRecHits);
  }

  // put the collection of recunstructed hits in the event
  LogInfo("EcalRecHitInfo") << "total # EB calibrated rechits: " << ebRecHits->size();
  LogInfo("EcalRecHitInfo") << "total # EE calibrated rechits: " << eeRecHits->size();

  evt.put(std::move(ebRecHits), ebRechitCollection_);
  evt.put(std::move(eeRecHits), eeRechitCollection_);
}

void EcalRecHitProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
  edm::ParameterSetDescription desc;
  desc.add<bool>("recoverEEVFE", false);
  desc.add<std::string>("EErechitCollection", "EcalRecHitsEE");
  desc.add<bool>("recoverEBIsolatedChannels", false);
  desc.add<bool>("recoverEBVFE", false);
  desc.add<bool>("laserCorrection", true);
  desc.add<double>("EBLaserMIN", 0.5);
  desc.add<bool>("killDeadChannels", true);
  {
    std::vector<int> temp1;
    temp1.reserve(3);
    temp1.push_back(14);
    temp1.push_back(78);
    temp1.push_back(142);
    desc.add<std::vector<int>>("dbStatusToBeExcludedEB", temp1);
  }
  desc.add<edm::InputTag>("EEuncalibRecHitCollection",
                          edm::InputTag("ecalMultiFitUncalibRecHit", "EcalUncalibRecHitsEE"));
  {
    std::vector<int> temp1;
    temp1.reserve(3);
    temp1.push_back(14);
    temp1.push_back(78);
    temp1.push_back(142);
    desc.add<std::vector<int>>("dbStatusToBeExcludedEE", temp1);
  }
  desc.add<double>("EELaserMIN", 0.5);
  desc.add<edm::InputTag>("ebFEToBeRecovered", edm::InputTag("ecalDetIdToBeRecovered", "ebFE"));
  {
    edm::ParameterSetDescription psd0;
    psd0.add<double>("e6e2thresh", 0.04);
    psd0.add<double>("tightenCrack_e6e2_double", 3);
    psd0.add<double>("e4e1Threshold_endcap", 0.3);
    psd0.add<double>("tightenCrack_e4e1_single", 3);
    psd0.add<double>("tightenCrack_e1_double", 2);
    psd0.add<double>("cThreshold_barrel", 4);
    psd0.add<double>("e4e1Threshold_barrel", 0.08);
    psd0.add<double>("tightenCrack_e1_single", 2);
    psd0.add<double>("e4e1_b_barrel", -0.024);
    psd0.add<double>("e4e1_a_barrel", 0.04);
    psd0.add<double>("ignoreOutOfTimeThresh", 1000000000.0);
    psd0.add<double>("cThreshold_endcap", 15);
    psd0.add<double>("e4e1_b_endcap", -0.0125);
    psd0.add<double>("e4e1_a_endcap", 0.02);
    psd0.add<double>("cThreshold_double", 10);
    desc.add<edm::ParameterSetDescription>("cleaningConfig", psd0);
  }
  desc.add<double>("logWarningEtThreshold_EE_FE", 50);
  desc.add<edm::InputTag>("eeDetIdToBeRecovered", edm::InputTag("ecalDetIdToBeRecovered", "eeDetId"));
  desc.add<bool>("recoverEBFE", true);
  desc.add<edm::InputTag>("eeFEToBeRecovered", edm::InputTag("ecalDetIdToBeRecovered", "eeFE"));
  desc.add<edm::InputTag>("ebDetIdToBeRecovered", edm::InputTag("ecalDetIdToBeRecovered", "ebDetId"));
  desc.add<double>("singleChannelRecoveryThreshold", 8);
  desc.add<double>("sum8ChannelRecoveryThreshold", 0.);
  desc.add<edm::FileInPath>("bdtWeightFileNoCracks",
                            edm::FileInPath("RecoLocalCalo/EcalDeadChannelRecoveryAlgos/data/BDTWeights/"
                                            "bdtgAllRH_8GT700MeV_noCracks_ZskimData2017_v1.xml"));
  desc.add<edm::FileInPath>("bdtWeightFileCracks",
                            edm::FileInPath("RecoLocalCalo/EcalDeadChannelRecoveryAlgos/data/BDTWeights/"
                                            "bdtgAllRH_8GT700MeV_onlyCracks_ZskimData2017_v1.xml"));
  {
    std::vector<std::string> temp1;
    temp1.reserve(9);
    temp1.push_back("kNoisy");
    temp1.push_back("kNNoisy");
    temp1.push_back("kFixedG6");
    temp1.push_back("kFixedG1");
    temp1.push_back("kFixedG0");
    temp1.push_back("kNonRespondingIsolated");
    temp1.push_back("kDeadVFE");
    temp1.push_back("kDeadFE");
    temp1.push_back("kNoDataNoTP");
    desc.add<std::vector<std::string>>("ChannelStatusToBeExcluded", temp1);
  }
  desc.add<std::string>("EBrechitCollection", "EcalRecHitsEB");
  desc.add<edm::InputTag>("triggerPrimitiveDigiCollection", edm::InputTag("ecalDigis", "EcalTriggerPrimitives"));
  desc.add<bool>("recoverEEFE", true);
  desc.add<std::string>("singleChannelRecoveryMethod", "NeuralNetworks");
  desc.add<double>("EBLaserMAX", 3.0);
  {
    edm::ParameterSetDescription psd0;
    {
      std::vector<std::string> temp2;
      temp2.reserve(4);
      temp2.push_back("kOk");
      temp2.push_back("kDAC");
      temp2.push_back("kNoLaser");
      temp2.push_back("kNoisy");
      psd0.add<std::vector<std::string>>("kGood", temp2);
    }
    {
      std::vector<std::string> temp2;
      temp2.reserve(3);
      temp2.push_back("kFixedG0");
      temp2.push_back("kNonRespondingIsolated");
      temp2.push_back("kDeadVFE");
      psd0.add<std::vector<std::string>>("kNeighboursRecovered", temp2);
    }
    {
      std::vector<std::string> temp2;
      temp2.reserve(1);
      temp2.push_back("kNoDataNoTP");
      psd0.add<std::vector<std::string>>("kDead", temp2);
    }
    {
      std::vector<std::string> temp2;
      temp2.reserve(3);
      temp2.push_back("kNNoisy");
      temp2.push_back("kFixedG6");
      temp2.push_back("kFixedG1");
      psd0.add<std::vector<std::string>>("kNoisy", temp2);
    }
    {
      std::vector<std::string> temp2;
      temp2.reserve(1);
      temp2.push_back("kDeadFE");
      psd0.add<std::vector<std::string>>("kTowerRecovered", temp2);
    }
    desc.add<edm::ParameterSetDescription>("flagsMapDBReco", psd0);
  }
  desc.add<edm::InputTag>("EBuncalibRecHitCollection",
                          edm::InputTag("ecalMultiFitUncalibRecHit", "EcalUncalibRecHitsEB"));
  desc.add<std::string>("algoRecover", "EcalRecHitWorkerRecover");
  desc.add<std::string>("algo", "EcalRecHitWorkerSimple");
  desc.add<double>("EELaserMAX", 8.0);
  desc.add<double>("logWarningEtThreshold_EB_FE", 50);
  desc.add<bool>("recoverEEIsolatedChannels", false);
  desc.add<bool>("skipTimeCalib", false);
  descriptions.add("ecalRecHit", desc);
}

#include "FWCore/Framework/interface/MakerMacros.h"
DEFINE_FWK_MODULE(EcalRecHitProducer);