EcalRecProducers/plugins/EcalTPSkimmer.cc

0001 /** \class EcalTPSkimmer
0002  *   produce a subset of TP information
0003  *
0004  *  \author Federico Ferri, CEA/Saclay Irfu/SPP
0005  *
0006  **/
0007
0008 #include "RecoLocalCalo/EcalRecProducers/plugins/EcalTPSkimmer.h"
0009
0010 #include "DataFormats/Common/interface/Handle.h"
0011 #include "FWCore/Framework/interface/ESHandle.h"
0012
0013 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0014
0015 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0016 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0017 #include "DataFormats/EcalDetId/interface/EcalTrigTowerDetId.h"
0018
0019 EcalTPSkimmer::EcalTPSkimmer(const edm::ParameterSet& ps) {
0020   skipModule_ = ps.getParameter<bool>("skipModule");
0021
0022   doBarrel_ = ps.getParameter<bool>("doBarrel");
0023   doEndcap_ = ps.getParameter<bool>("doEndcap");
0024
0025   chStatusToSelectTP_ = ps.getParameter<std::vector<uint32_t> >("chStatusToSelectTP");
0026
0027   tpOutputCollection_ = ps.getParameter<std::string>("tpOutputCollection");
0028   tpInputToken_ = consumes<EcalTrigPrimDigiCollection>(ps.getParameter<edm::InputTag>("tpInputCollection"));
0029   ttMapToken_ = esConsumes<EcalTrigTowerConstituentsMap, IdealGeometryRecord>();
0030   if (not skipModule_) {
0031     chStatusToken_ = esConsumes<EcalChannelStatus, EcalChannelStatusRcd>();
0032   }
0033   produces<EcalTrigPrimDigiCollection>(tpOutputCollection_);
0034 }
0035
0036 EcalTPSkimmer::~EcalTPSkimmer() {}
0037
0038 void EcalTPSkimmer::produce(edm::Event& evt, const edm::EventSetup& es) {
0039   insertedTP_.clear();
0040
0041   using namespace edm;
0042
0043   ttMap_ = es.getHandle(ttMapToken_);
0044
0045   // collection of rechits to put in the event
0046   auto tpOut = std::make_unique<EcalTrigPrimDigiCollection>();
0047
0048   if (skipModule_) {
0049     evt.put(std::move(tpOut), tpOutputCollection_);
0050     return;
0051   }
0052
0053   edm::ESHandle<EcalChannelStatus> chStatus = es.getHandle(chStatusToken_);
0054
0055   edm::Handle<EcalTrigPrimDigiCollection> tpIn;
0056   evt.getByToken(tpInputToken_, tpIn);
0057
0058   if (doBarrel_) {
0059     EcalChannelStatusMap::const_iterator chit;
0060     uint16_t code = 0;
0061     for (int i = 0; i < EBDetId::kSizeForDenseIndexing; ++i) {
0062       if (!EBDetId::validDenseIndex(i))
0063         continue;
0064       EBDetId id = EBDetId::detIdFromDenseIndex(i);
0065       chit = chStatus->find(id);
0066       // check if the channel status means TP to be kept
0067       if (chit != chStatus->end()) {
0068         code = (*chit).getStatusCode();
0069         if (std::find(chStatusToSelectTP_.begin(), chStatusToSelectTP_.end(), code) != chStatusToSelectTP_.end()) {
0070           // retrieve the TP DetId
0071           EcalTrigTowerDetId ttDetId(((EBDetId)id).tower());
0072           // insert the TP if not done already
0073           if (!alreadyInserted(ttDetId))
0074             insertTP(ttDetId, tpIn, *tpOut);
0075         }
0076       } else {
0077         edm::LogError("EcalDetIdToBeRecoveredProducer") << "No channel status found for xtal " << id.rawId()
0078                                                         << "! something wrong with EcalChannelStatus in your DB? ";
0079       }
0080     }
0081   }
0082
0083   if (doEndcap_) {
0084     EcalChannelStatusMap::const_iterator chit;
0085     uint16_t code = 0;
0086     for (int i = 0; i < EEDetId::kSizeForDenseIndexing; ++i) {
0087       if (!EEDetId::validDenseIndex(i))
0088         continue;
0089       EEDetId id = EEDetId::detIdFromDenseIndex(i);
0090       chit = chStatus->find(id);
0091       // check if the channel status means TP to be kept
0092       if (chit != chStatus->end()) {
0093         code = (*chit).getStatusCode();
0094         if (std::find(chStatusToSelectTP_.begin(), chStatusToSelectTP_.end(), code) != chStatusToSelectTP_.end()) {
0095           // retrieve the TP DetId
0096           EcalTrigTowerDetId ttDetId = ttMap_->towerOf(id);
0097           // insert the TP if not done already
0098           if (!alreadyInserted(ttDetId))
0099             insertTP(ttDetId, tpIn, *tpOut);
0100         }
0101       } else {
0102         edm::LogError("EcalDetIdToBeRecoveredProducer") << "No channel status found for xtal " << id.rawId()
0103                                                         << "! something wrong with EcalChannelStatus in your DB? ";
0104       }
0105     }
0106   }
0107
0108   // put the collection of reconstructed hits in the event
0109   LogInfo("EcalTPSkimmer") << "total # of TP inserted: " << tpOut->size();
0110
0111   evt.put(std::move(tpOut), tpOutputCollection_);
0112 }
0113
0114 bool EcalTPSkimmer::alreadyInserted(EcalTrigTowerDetId ttId) { return (insertedTP_.find(ttId) != insertedTP_.end()); }
0115
0116 void EcalTPSkimmer::insertTP(EcalTrigTowerDetId ttId,
0117                              edm::Handle<EcalTrigPrimDigiCollection>& tpIn,
0118                              EcalTrigPrimDigiCollection& tpOut) {
0119   EcalTrigPrimDigiCollection::const_iterator tpIt = tpIn->find(ttId);
0120   if (tpIt != tpIn->end()) {
0121     tpOut.push_back(*tpIt);
0122     insertedTP_.insert(ttId);
0123   }
0124 }
0125
0126 #include "FWCore/Framework/interface/MakerMacros.h"
0127 DEFINE_FWK_MODULE(EcalTPSkimmer);