HcalAlCaRecoProducers/plugins/AlCaHcalNoiseProducer.cc

0001 /*
0002 Original author Grigory Safronov
0003
0004 27/03/09 - compilation from :
0005 HLTrigger/special/src/HLTHcalNoiseFilter.cc
0006 Calibration/HcalAlCaRecoProducers/src/AlCaEcalHcalReadoutsProducer.cc
0007 Calibration/HcalIsolatedTrackReco/src/SubdetFEDSelector.cc
0008
0009 */
0010
0011 // -*- C++ -*-
0012
0013 // system include files
0014 #include <memory>
0015 #include <string>
0016 // user include files
0017
0018 #include "FWCore/Framework/interface/Frameworkfwd.h"
0019 #include "FWCore/Framework/interface/one/EDProducer.h"
0020 #include "FWCore/Framework/interface/Event.h"
0021 #include "FWCore/Framework/interface/EventSetup.h"
0022 #include "FWCore/Framework/interface/MakerMacros.h"
0023 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0024 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0025 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0026
0027 #include "DataFormats/Common/interface/Ref.h"
0028 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
0029 #include "DataFormats/DetId/interface/DetId.h"
0030
0031 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0032 #include "DataFormats/CaloTowers/interface/CaloTowerDetId.h"
0033 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0034 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0035
0036 #include "DataFormats/CaloTowers/interface/CaloTower.h"
0037 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
0038 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0039 #include "DataFormats/JetReco/interface/CaloJet.h"
0040 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0041 #include "DataFormats/Math/interface/deltaR.h"
0042 #include "DataFormats/METReco/interface/CaloMET.h"
0043 #include "DataFormats/METReco/interface/CaloMETCollection.h"
0044
0045 #include "DataFormats/FEDRawData/interface/FEDRawData.h"
0046 #include "DataFormats/FEDRawData/interface/FEDNumbering.h"
0047 #include "DataFormats/FEDRawData/interface/FEDRawDataCollection.h"
0048
0049 #include "EventFilter/RawDataCollector/interface/RawDataFEDSelector.h"
0050 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0051 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0052
0053 //
0054 // class decleration
0055 //
0056
0057 class AlCaHcalNoiseProducer : public edm::one::EDProducer<> {
0058 public:
0059   explicit AlCaHcalNoiseProducer(const edm::ParameterSet&);
0060   ~AlCaHcalNoiseProducer() override = default;
0061
0062   void produce(edm::Event&, const edm::EventSetup&) override;
0063   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0064
0065 private:
0066   // ----------member data ---------------------------
0067
0068   bool useMet_;
0069   bool useJet_;
0070   double MetCut_;
0071   double JetMinE_;
0072   double JetHCALminEnergyFraction_;
0073   int nAnomalousEvents;
0074   int nEvents;
0075
0076   std::vector<edm::InputTag> ecalLabels_;
0077
0078   edm::EDGetTokenT<reco::CaloJetCollection> tok_jets_;
0079   edm::EDGetTokenT<reco::CaloMETCollection> tok_met_;
0080   edm::EDGetTokenT<CaloTowerCollection> tok_tower_;
0081
0082   edm::EDGetTokenT<HBHERecHitCollection> tok_hbhe_;
0083   edm::EDGetTokenT<HORecHitCollection> tok_ho_;
0084   edm::EDGetTokenT<HFRecHitCollection> tok_hf_;
0085
0086   edm::EDGetTokenT<EcalRecHitCollection> tok_ps_;
0087   edm::EDGetTokenT<FEDRawDataCollection> tok_raw_;
0088   std::vector<edm::EDGetTokenT<EcalRecHitCollection>> toks_ecal_;
0089 };
0090
0091 AlCaHcalNoiseProducer::AlCaHcalNoiseProducer(const edm::ParameterSet& iConfig) {
0092   tok_jets_ = consumes<reco::CaloJetCollection>(iConfig.getParameter<edm::InputTag>("JetSource"));
0093   tok_met_ = consumes<reco::CaloMETCollection>(iConfig.getParameter<edm::InputTag>("MetSource"));
0094   tok_tower_ = consumes<CaloTowerCollection>(iConfig.getParameter<edm::InputTag>("TowerSource"));
0095   useMet_ = iConfig.getParameter<bool>("UseMET");
0096   useJet_ = iConfig.getParameter<bool>("UseJet");
0097   MetCut_ = iConfig.getParameter<double>("MetCut");
0098   JetMinE_ = iConfig.getParameter<double>("JetMinE");
0099   JetHCALminEnergyFraction_ = iConfig.getParameter<double>("JetHCALminEnergyFraction");
0100
0101   tok_ho_ = consumes<HORecHitCollection>(iConfig.getParameter<edm::InputTag>("hoInput"));
0102   tok_hf_ = consumes<HFRecHitCollection>(iConfig.getParameter<edm::InputTag>("hfInput"));
0103   tok_hbhe_ = consumes<HBHERecHitCollection>(iConfig.getParameter<edm::InputTag>("hbheInput"));
0104   ecalLabels_ = iConfig.getParameter<std::vector<edm::InputTag>>("ecalInputs");
0105   tok_ps_ = consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("ecalPSInput"));
0106   tok_raw_ = consumes<FEDRawDataCollection>(iConfig.getParameter<edm::InputTag>("rawInput"));
0107
0108   const unsigned nLabels = ecalLabels_.size();
0109   for (unsigned i = 0; i != nLabels; i++)
0110     toks_ecal_.push_back(consumes<EcalRecHitCollection>(ecalLabels_[i]));
0111
0112   //register products
0113   produces<HBHERecHitCollection>("HBHERecHitCollectionFHN");
0114   produces<HORecHitCollection>("HORecHitCollectionFHN");
0115   produces<HFRecHitCollection>("HFRecHitCollectionFHN");
0116
0117   produces<EcalRecHitCollection>("EcalRecHitCollectionFHN");
0118   produces<EcalRecHitCollection>("PSEcalRecHitCollectionFHN");
0119
0120   produces<FEDRawDataCollection>("HcalFEDsFHN");
0121 }
0122
0123 // ------------ method called to produce the data  ------------
0124 void AlCaHcalNoiseProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0125   bool acceptEvent = false;
0126
0127   // filtering basing on HLTrigger/special/src/HLTHcalNoiseFilter.cc:
0128
0129   bool isAnomalous_BasedOnMET = false;
0130   bool isAnomalous_BasedOnEnergyFraction = false;
0131
0132   if (useMet_) {
0133     edm::Handle<reco::CaloMETCollection> metHandle = iEvent.getHandle(tok_met_);
0134     const reco::CaloMETCollection* metCol = metHandle.product();
0135     const reco::CaloMET met = metCol->front();
0136
0137     if (met.pt() > MetCut_)
0138       isAnomalous_BasedOnMET = true;
0139   }
0140
0141   if (useJet_) {
0142     edm::Handle<reco::CaloJetCollection> calojetHandle = iEvent.getHandle(tok_jets_);
0143     edm::Handle<CaloTowerCollection> towerHandle = iEvent.getHandle(tok_tower_);
0144
0145     std::vector<CaloTower> TowerContainer;
0146     std::vector<reco::CaloJet> JetContainer;
0147     TowerContainer.clear();
0148     JetContainer.clear();
0149     CaloTower seedTower;
0150     nEvents++;
0151     for (reco::CaloJetCollection::const_iterator calojetIter = calojetHandle->begin();
0152          calojetIter != calojetHandle->end();
0153          ++calojetIter) {
0154       if (((calojetIter->et()) * cosh(calojetIter->eta()) > JetMinE_) &&
0155           (calojetIter->energyFractionHadronic() > JetHCALminEnergyFraction_)) {
0156         JetContainer.push_back(*calojetIter);
0157         double maxTowerE = 0.0;
0158         for (CaloTowerCollection::const_iterator kal = towerHandle->begin(); kal != towerHandle->end(); kal++) {
0159           double dR = deltaR((*calojetIter).eta(), (*calojetIter).phi(), (*kal).eta(), (*kal).phi());
0160           if ((dR < 0.50) && (kal->p() > maxTowerE)) {
0161             maxTowerE = kal->p();
0162             seedTower = *kal;
0163           }
0164         }
0165         TowerContainer.push_back(seedTower);
0166       }
0167     }
0168     if (!JetContainer.empty()) {
0169       nAnomalousEvents++;
0170       isAnomalous_BasedOnEnergyFraction = true;
0171     }
0172   }
0173
0174   acceptEvent = ((useMet_ && isAnomalous_BasedOnMET) || (useJet_ && isAnomalous_BasedOnEnergyFraction));
0175
0176   ////////////////
0177
0178   //Create empty output collections
0179
0180   auto miniHBHERecHitCollection = std::make_unique<HBHERecHitCollection>();
0181   auto miniHORecHitCollection = std::make_unique<HORecHitCollection>();
0182   auto miniHFRecHitCollection = std::make_unique<HFRecHitCollection>();
0183
0184   auto outputEColl = std::make_unique<EcalRecHitCollection>();
0185   auto outputESColl = std::make_unique<EcalRecHitCollection>();
0186
0187   auto outputFEDs = std::make_unique<FEDRawDataCollection>();
0188
0189   // if good event get and save all colletions
0190   if (acceptEvent) {
0191     edm::Handle<HBHERecHitCollection> hbhe = iEvent.getHandle(tok_hbhe_);
0192     edm::Handle<HORecHitCollection> ho = iEvent.getHandle(tok_ho_);
0193     edm::Handle<HFRecHitCollection> hf = iEvent.getHandle(tok_hf_);
0194
0195     edm::Handle<EcalRecHitCollection> pRecHits = iEvent.getHandle(tok_ps_);
0196
0197     // temporary collection of EB+EE recHits
0198
0199     auto tmpEcalRecHitCollection = std::make_unique<EcalRecHitCollection>();
0200
0201     std::vector<edm::EDGetTokenT<EcalRecHitCollection>>::const_iterator i;
0202     for (i = toks_ecal_.begin(); i != toks_ecal_.end(); i++) {
0203       edm::Handle<EcalRecHitCollection> ec = iEvent.getHandle(*i);
0204       for (EcalRecHitCollection::const_iterator recHit = (*ec).begin(); recHit != (*ec).end(); ++recHit) {
0205         tmpEcalRecHitCollection->push_back(*recHit);
0206       }
0207     }
0208
0209     //////////
0210
0211     //////// write HCAL collections:
0212     const HBHERecHitCollection Hithbhe = *(hbhe.product());
0213     for (HBHERecHitCollection::const_iterator hbheItr = Hithbhe.begin(); hbheItr != Hithbhe.end(); hbheItr++) {
0214       miniHBHERecHitCollection->push_back(*hbheItr);
0215     }
0216     const HORecHitCollection Hitho = *(ho.product());
0217     for (HORecHitCollection::const_iterator hoItr = Hitho.begin(); hoItr != Hitho.end(); hoItr++) {
0218       miniHORecHitCollection->push_back(*hoItr);
0219     }
0220
0221     const HFRecHitCollection Hithf = *(hf.product());
0222     for (HFRecHitCollection::const_iterator hfItr = Hithf.begin(); hfItr != Hithf.end(); hfItr++) {
0223       miniHFRecHitCollection->push_back(*hfItr);
0224     }
0225     /////
0226
0227     ///// write ECAL
0228     for (std::vector<EcalRecHit>::const_iterator ehit = tmpEcalRecHitCollection->begin();
0229          ehit != tmpEcalRecHitCollection->end();
0230          ehit++) {
0231       outputEColl->push_back(*ehit);
0232     }
0233     /////////
0234
0235     // write PS
0236     const EcalRecHitCollection& psrechits = *(pRecHits.product());
0237
0238     for (EcalRecHitCollection::const_iterator i = psrechits.begin(); i != psrechits.end(); i++) {
0239       outputESColl->push_back(*i);
0240     }
0241
0242     // get HCAL FEDs
0243     edm::Handle<FEDRawDataCollection> rawIn = iEvent.getHandle(tok_raw_);
0244
0245     std::vector<int> selFEDs;
0246     for (int i = FEDNumbering::MINHCALFEDID; i <= FEDNumbering::MAXHCALFEDID; i++) {
0247       selFEDs.push_back(i);
0248     }
0249     ////////////
0250
0251     // Copying FEDs :
0252     const FEDRawDataCollection* rdc = rawIn.product();
0253
0254     //   if ( ( rawData[i].provenance()->processName() != e.processHistory().rbegin()->processName() ) )
0255     //       continue ; // skip all raw collections not produced by the current process
0256
0257     for (int j = 0; j <= FEDNumbering::MAXFEDID; ++j) {
0258       bool rightFED = false;
0259       for (uint32_t k = 0; k < selFEDs.size(); k++) {
0260         if (j == selFEDs[k]) {
0261           rightFED = true;
0262         }
0263       }
0264       if (!rightFED)
0265         continue;
0266       const FEDRawData& fedData = rdc->FEDData(j);
0267       size_t size = fedData.size();
0268
0269       if (size > 0) {
0270         // this fed has data -- lets copy it
0271         FEDRawData& fedDataProd = outputFEDs->FEDData(j);
0272         if (fedDataProd.size() != 0) {
0273           edm::LogWarning("HcalNoise") << " More than one FEDRawDataCollection with data in FED " << j
0274                                        << " Skipping the 2nd\n";
0275           continue;
0276         }
0277         fedDataProd.resize(size);
0278         unsigned char* dataProd = fedDataProd.data();
0279         const unsigned char* data = fedData.data();
0280         for (unsigned int k = 0; k < size; ++k) {
0281           dataProd[k] = data[k];
0282         }
0283       }
0284     }
0285     //////////////////////
0286   }
0287
0288   //Put selected information in the event
0289   iEvent.put(std::move(miniHBHERecHitCollection), "HBHERecHitCollectionFHN");
0290   iEvent.put(std::move(miniHORecHitCollection), "HORecHitCollectionFHN");
0291   iEvent.put(std::move(miniHFRecHitCollection), "HFRecHitCollectionFHN");
0292   iEvent.put(std::move(outputEColl), "EcalRecHitCollectionFHN");
0293   iEvent.put(std::move(outputESColl), "PSEcalRecHitCollectionFHN");
0294   iEvent.put(std::move(outputFEDs), "HcalFEDsFHN");
0295 }
0296
0297 void AlCaHcalNoiseProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0298   edm::ParameterSetDescription desc;
0299   desc.add<edm::InputTag>("JetSource", edm::InputTag("iterativeCone5CaloJets"));
0300   desc.add<edm::InputTag>("MetSource", edm::InputTag("met"));
0301   desc.add<edm::InputTag>("TowerSource", edm::InputTag("towerMaker"));
0302   desc.add<bool>("UseJet", true);
0303   desc.add<bool>("UseMET", false);
0304   desc.add<double>("MetCut", 0);
0305   desc.add<double>("JetMinE", 20);
0306   desc.add<double>("JetHCALminEnergyFraction", 0.98);
0307   desc.add<edm::InputTag>("hbheInput", edm::InputTag("hbhereco"));
0308   desc.add<edm::InputTag>("hfInput", edm::InputTag("hfreco"));
0309   desc.add<edm::InputTag>("hoInput", edm::InputTag("horeco"));
0310   std::vector<edm::InputTag> inputs = {edm::InputTag("ecalRecHit", "EcalRecHitsEB"),
0311                                        edm::InputTag("ecalRecHit", "EcalRecHitsEE")};
0312   desc.add<std::vector<edm::InputTag>>("ecalInputs", inputs);
0313   desc.add<edm::InputTag>("ecalPSInput", edm::InputTag("ecalPreshowerRecHit", "EcalRecHitsES"));
0314   desc.add<edm::InputTag>("rawInput", edm::InputTag("rawDataCollector"));
0315   descriptions.add("alcaHcalNoiseProducer", desc);
0316 }
0317
0318 #include "FWCore/Framework/interface/MakerMacros.h"
0319
0320 DEFINE_FWK_MODULE(AlCaHcalNoiseProducer);