HiEgammaAlgos/plugins/HiSpikeCleaner.cc

0001 // -*- C++ -*-
0002 //111
0003 // Package:    HiSpikeCleaner
0004 // Class:      HiSpikeCleaner
0005 //
0006 /**\class HiSpikeCleaner HiSpikeCleaner.cc RecoHI/HiSpikeCleaner/src/HiSpikeCleaner.cc
0007
0008  Description: [one line class summary]
0009
0010  Implementation:
0011      [Notes on implementation]
0012 */
0013 //
0014 // Original Author:  Yong Kim,32 4-A08,+41227673039,
0015 //         Created:  Mon Nov  1 18:22:21 CET 2010
0016 //
0017 //
0018
0019 // system include files
0020 #include <memory>
0021
0022 // user include files
0023 #include "FWCore/Framework/interface/Frameworkfwd.h"
0024 #include "FWCore/Framework/interface/stream/EDProducer.h"
0025 #include "FWCore/Framework/interface/EventSetup.h"
0026
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/Framework/interface/MakerMacros.h"
0029
0030 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0031 #include "FWCore/Utilities/interface/InputTag.h"
0032
0033 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0034 #include "DataFormats/EgammaCandidates/interface/PhotonFwd.h"
0035
0036 #include "RecoEcal/EgammaCoreTools/interface/EcalTools.h"
0037 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgoRcd.h"
0038 #include "RecoLocalCalo/EcalRecAlgos/interface/EcalSeverityLevelAlgo.h"
0039 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h"
0040
0041 //
0042 // class declaration
0043 //
0044
0045 class HiSpikeCleaner : public edm::stream::EDProducer<> {
0046 public:
0047   explicit HiSpikeCleaner(const edm::ParameterSet&);
0048   ~HiSpikeCleaner() override;
0049
0050 private:
0051   void produce(edm::Event&, const edm::EventSetup&) override;
0052
0053   // ----------member data ---------------------------
0054
0055   edm::EDGetTokenT<reco::SuperClusterCollection> sCInputProducerToken_;
0056   edm::EDGetTokenT<EcalRecHitCollection> rHInputProducerBToken_;
0057   edm::EDGetTokenT<EcalRecHitCollection> rHInputProducerEToken_;
0058   edm::ESGetToken<EcalSeverityLevelAlgo, EcalSeverityLevelAlgoRcd> ecalSevLvlAlgoToken_;
0059   const EcalClusterLazyTools::ESGetTokens ecalClusterToolsESGetTokens_;
0060
0061   std::string outputCollection_;
0062   double TimingCut_;
0063   double swissCutThr_;
0064   double etCut_;
0065 };
0066
0067 HiSpikeCleaner::HiSpikeCleaner(const edm::ParameterSet& iConfig)
0068     : ecalSevLvlAlgoToken_(esConsumes()), ecalClusterToolsESGetTokens_{consumesCollector()} {
0069   //register your products
0070   /* Examples
0071    produces<ExampleData2>();
0072
0073    //if do put with a label
0074    produces<ExampleData2>("label");
0075 */
0076   //now do what ever other initialization is needed
0077
0078   rHInputProducerBToken_ = consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("recHitProducerBarrel"));
0079   rHInputProducerEToken_ = consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("recHitProducerEndcap"));
0080
0081   sCInputProducerToken_ =
0082       consumes<reco::SuperClusterCollection>(iConfig.getParameter<edm::InputTag>("originalSuperClusterProducer"));
0083   TimingCut_ = iConfig.getUntrackedParameter<double>("TimingCut", 4.0);
0084   swissCutThr_ = iConfig.getUntrackedParameter<double>("swissCutThr", 0.95);
0085   etCut_ = iConfig.getParameter<double>("etCut");
0086
0087   outputCollection_ = iConfig.getParameter<std::string>("outputColl");
0088   produces<reco::SuperClusterCollection>(outputCollection_);
0089 }
0090
0091 HiSpikeCleaner::~HiSpikeCleaner() {
0092   // do anything here that needs to be done at desctruction time
0093   // (e.g. close files, deallocate resources etc.)
0094 }
0095
0096 //
0097 // member functions
0098 //
0099
0100 // ------------ method called to produce the data  ------------
0101 void HiSpikeCleaner::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0102   using namespace edm;
0103
0104   // Get raw SuperClusters from the event
0105   Handle<reco::SuperClusterCollection> pRawSuperClusters;
0106   iEvent.getByToken(sCInputProducerToken_, pRawSuperClusters);
0107
0108   // Get the RecHits from the event
0109   Handle<EcalRecHitCollection> pRecHitsB;
0110   iEvent.getByToken(rHInputProducerBToken_, pRecHitsB);
0111
0112   // Get the RecHits from the event
0113   Handle<EcalRecHitCollection> pRecHitsE;
0114   iEvent.getByToken(rHInputProducerEToken_, pRecHitsE);
0115
0116   // get the channel status from the DB
0117   //   edm::ESHandle<EcalChannelStatus> chStatus;
0118   //   iSetup.get<EcalChannelStatusRcd>().get(chStatus);
0119
0120   auto const& ecalSevLvlAlgo = iSetup.getData(ecalSevLvlAlgoToken_);
0121
0122   // Create a pointer to the RecHits and raw SuperClusters
0123   const reco::SuperClusterCollection* rawClusters = pRawSuperClusters.product();
0124
0125   EcalClusterLazyTools lazyTool(
0126       iEvent, ecalClusterToolsESGetTokens_.get(iSetup), rHInputProducerBToken_, rHInputProducerEToken_);
0127
0128   // Define a collection of corrected SuperClusters to put back into the event
0129   auto corrClusters = std::make_unique<reco::SuperClusterCollection>();
0130
0131   //  Loop over raw clusters and make corrected ones
0132   reco::SuperClusterCollection::const_iterator aClus;
0133   for (aClus = rawClusters->begin(); aClus != rawClusters->end(); aClus++) {
0134     double theEt = aClus->energy() / cosh(aClus->eta());
0135     //   std::cout << " et of SC = " << theEt << std::endl;
0136
0137     if (theEt < etCut_)
0138       continue;  // cut off low pT superclusters
0139
0140     bool flagS = true;
0141     float swissCrx(0);
0142
0143     const reco::CaloClusterPtr seed = aClus->seed();
0144     DetId id = lazyTool.getMaximum(*seed).first;
0145     const EcalRecHitCollection& rechits = *pRecHitsB;
0146     EcalRecHitCollection::const_iterator it = rechits.find(id);
0147
0148     if (it != rechits.end()) {
0149       ecalSevLvlAlgo.severityLevel(id, rechits);
0150       swissCrx = EcalTools::swissCross(id, rechits, 0., true);
0151       //        std::cout << "swissCross = " << swissCrx <<std::endl;
0152       // std::cout << " timing = " << it->time() << std::endl;
0153     }
0154
0155     if (fabs(it->time()) > TimingCut_) {
0156       flagS = false;
0157       //        std::cout << " timing = " << it->time() << std::endl;
0158       //   std::cout << " timing is bad........" << std::endl;
0159     }
0160     if (swissCrx > (float)swissCutThr_) {
0161       flagS = false;  // swissCross cut
0162                       //        std::cout << "swissCross = " << swissCrx <<std::endl;
0163                       //   std::cout << " removed by swiss cross cut" << std::endl;
0164     }
0165     // - kGood        --> good channel
0166     // - kProblematic --> problematic (e.g. noisy)
0167     // - kRecovered   --> recovered (e.g. an originally dead or saturated)
0168     // - kTime        --> the channel is out of time (e.g. spike)
0169     // - kWeird       --> weird (e.g. spike)
0170     // - kBad         --> bad, not suitable to be used in the reconstruction
0171     //   enum EcalSeverityLevel { kGood=0, kProblematic, kRecovered, kTime, kWeird, kBad };
0172
0173     reco::SuperCluster newClus;
0174     if (flagS == true)
0175       newClus = *aClus;
0176     else
0177       continue;
0178     corrClusters->push_back(newClus);
0179   }
0180
0181   // Put collection of corrected SuperClusters into the event
0182   iEvent.put(std::move(corrClusters), outputCollection_);
0183 }
0184
0185 //define this as a plug-in
0186 DEFINE_FWK_MODULE(HiSpikeCleaner);