special/plugins/EcalMIPRecHitFilter.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    EcalMIPRecHitFilter
0004 // Class:      EcalMIPRecHitFilter
0005 //
0006 /**\class EcalMIPRecHitFilter EcalMIPRecHitFilter.cc Work/EcalMIPRecHitFilter/src/EcalMIPRecHitFilter.cc
0007
0008  Description: <one line class summary>
0009
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Giovanni FRANZONI
0015 //         Created:  Wed Sep 19 16:21:29 CEST 2007
0016 //
0017 //
0018
0019 // system include files
0020 #include <memory>
0021 #include <vector>
0022
0023 // user include files
0024 #include "FWCore/Framework/interface/Frameworkfwd.h"
0025 #include "FWCore/Framework/interface/global/EDFilter.h"
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/ESHandle.h"
0028
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0031 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0032
0033 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0034
0035 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
0036 #include "DataFormats/EcalRecHit/interface/EcalUncalibratedRecHit.h"
0037 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0038
0039 #include "DataFormats/DetId/interface/DetId.h"
0040 #include "Geometry/CaloTopology/interface/CaloTopology.h"
0041 #include "Geometry/Records/interface/CaloTopologyRecord.h"
0042
0043 #include "CondFormats/EcalObjects/interface/EcalIntercalibConstants.h"
0044 #include "CondFormats/DataRecord/interface/EcalIntercalibConstantsRcd.h"
0045 #include "CalibCalorimetry/EcalLaserCorrection/interface/EcalLaserDbService.h"
0046 #include "CalibCalorimetry/EcalLaserCorrection/interface/EcalLaserDbRecord.h"
0047
0048 #include "CondFormats/EcalObjects/interface/EcalADCToGeVConstant.h"
0049 #include "CondFormats/DataRecord/interface/EcalADCToGeVConstantRcd.h"
0050
0051 //
0052 // class declaration
0053 //
0054
0055 class EcalMIPRecHitFilter : public edm::global::EDFilter<> {
0056 public:
0057   explicit EcalMIPRecHitFilter(const edm::ParameterSet&);
0058
0059   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0060
0061 private:
0062   bool filter(edm::StreamID, edm::Event&, edm::EventSetup const&) const override;
0063
0064   // ----------member data ---------------------------
0065   const edm::ESGetToken<CaloTopology, CaloTopologyRecord> caloTopologyRecordToken_;
0066   const edm::ESGetToken<EcalIntercalibConstants, EcalIntercalibConstantsRcd> ecalIntercalibConstantsRcdToken_;
0067   const edm::ESGetToken<EcalLaserDbService, EcalLaserDbRecord> ecalLaserDbRecordToken_;
0068   const edm::ESGetToken<EcalADCToGeVConstant, EcalADCToGeVConstantRcd> ecalADCToGeVConstantRcdToken_;
0069
0070   const edm::EDGetTokenT<EcalRecHitCollection> EcalRecHitToken_;
0071   const double minAmp1_;
0072   const double minAmp2_;
0073   const double minSingleAmp_;
0074   const std::vector<int> maskedList_;
0075   const int side_;
0076 };
0077
0078 //
0079 // constructors and destructor
0080 //
0081 EcalMIPRecHitFilter::EcalMIPRecHitFilter(const edm::ParameterSet& iConfig)
0082     : caloTopologyRecordToken_(esConsumes()),
0083       ecalIntercalibConstantsRcdToken_(esConsumes()),
0084       ecalLaserDbRecordToken_(esConsumes()),
0085       ecalADCToGeVConstantRcdToken_(esConsumes()),
0086       EcalRecHitToken_(consumes<EcalRecHitCollection>(iConfig.getParameter<edm::InputTag>("EcalRecHitCollection"))),
0087       minAmp1_(iConfig.getUntrackedParameter<double>("AmpMinSeed", 0.063)),
0088       minAmp2_(iConfig.getUntrackedParameter<double>("AmpMin2", 0.045)),
0089       minSingleAmp_(iConfig.getUntrackedParameter<double>("SingleAmpMin", 0.108)),
0090       maskedList_(iConfig.getUntrackedParameter<std::vector<int>>(
0091           "maskedChannels", std::vector<int>{})),  // this is using the ashed index
0092       side_(iConfig.getUntrackedParameter<int>("side", 3)) {
0093   // now do what ever initialization is needed
0094 }
0095
0096 //
0097 // member functions
0098 //
0099
0100 // ------------ method called on each new Event  ------------
0101 bool EcalMIPRecHitFilter::filter(edm::StreamID, edm::Event& iEvent, edm::EventSetup const& iSetup) const {
0102   using namespace edm;
0103
0104   // getting very basic uncalRH
0105   Handle<EcalRecHitCollection> recHits;
0106   if (not iEvent.getByToken(EcalRecHitToken_, recHits)) {
0107     edm::EDConsumerBase::Labels labels;
0108     labelsForToken(EcalRecHitToken_, labels);
0109     LogWarning("EcalMIPRecHitFilter") << "InputTag:  label = \"" << labels.module << "\", instance = \""
0110                                       << labels.productInstance << "\", process = \"" << labels.process
0111                                       << "\" is not available";
0112     return false;
0113   }
0114
0115   auto const& caloTopo = iSetup.getHandle(caloTopologyRecordToken_);
0116
0117   // Intercalib constants
0118   auto const& ical = iSetup.getData(ecalIntercalibConstantsRcdToken_);
0119   const EcalIntercalibConstantMap& icalMap = ical.getMap();
0120
0121   auto const& pLaser = iSetup.getHandle(ecalLaserDbRecordToken_);
0122
0123   auto const& agc = iSetup.getData(ecalADCToGeVConstantRcdToken_);
0124   //std::cout << "Global EB ADC->GeV scale: " << agc.getEBValue() << " GeV/ADC count" ;
0125   float adcconst = agc.getEBValue();
0126
0127   bool thereIsSignal = false;
0128   // loop on  rechits
0129   for (auto hitItr = recHits->begin(); hitItr != recHits->end(); ++hitItr) {
0130     EcalRecHit const& hit = *hitItr;
0131
0132     // masking noisy channels //KEEP this for now, just in case a few show up
0133     auto result = std::find(maskedList_.begin(), maskedList_.end(), EBDetId(hit.id()).hashedIndex());
0134     if (result != maskedList_.end())
0135       // LogWarning("EcalFilter") << "skipping uncalRecHit for channel: " << ic << " with amplitude " << ampli_ ;
0136       continue;
0137
0138     float ampli_ = hit.energy();
0139     EBDetId ebDet = hit.id();
0140
0141     // find intercalib constant for this xtal
0142     auto icalit = icalMap.find(ebDet);
0143     EcalIntercalibConstant icalconst = 1.;
0144     if (icalit != icalMap.end()) {
0145       icalconst = (*icalit);
0146       //LogDebug("EcalRecHitDebug") << "Found intercalib for xtal " << EBDetId(it->id()).ic() << " " << icalconst ;
0147     } else {
0148       //edm::LogError("EcalRecHitError") << "No intercalib const found for xtal " << EBDetId(ebDet) << "! something wrong with EcalIntercalibConstants in your DB? " ;
0149     }
0150     float lasercalib = pLaser->getLaserCorrection(EBDetId(ebDet), iEvent.time());
0151
0152     ampli_ /=
0153         (icalconst * lasercalib * adcconst);  ///LASER and CALIB constants from the DB //PUT THRESHOLDS IN ADC... AGAIN.
0154     // seeking channels with signal and displaced jitter
0155     if (ampli_ >= minSingleAmp_) {
0156       //std::cout << " THIS AMPLITUDE WORKS " << ampli_ << std::endl;
0157       thereIsSignal = true;
0158       // LogWarning("EcalFilter")  << "at evet: " << iEvent.id().event()
0159       //                       << " and run: " << iEvent.id().run()
0160       //                       << " there is OUT OF TIME signal at chanel: " << ic
0161       //                       << " with amplitude " << ampli_  << " and max at: " << jitter_;
0162       break;
0163     }
0164
0165     //Check for more robust selection other than just single crystal cosmics
0166     if (ampli_ >= minAmp1_) {
0167       //std::cout << " THIS AMPLITUDE WORKS " << ampli_ << std::endl;
0168       std::vector<DetId> neighbors = caloTopo->getWindow(ebDet, side_, side_);
0169       float secondMin = 0.;
0170       for (std::vector<DetId>::const_iterator detitr = neighbors.begin(); detitr != neighbors.end(); ++detitr) {
0171         auto thishit = recHits->find((*detitr));
0172         if (thishit == recHits->end()) {
0173           //LogWarning("EcalMIPRecHitFilter") << "No RecHit available, for "<< EBDetId(*detitr);
0174           continue;
0175         }
0176         if ((*thishit).id() != ebDet) {
0177           float thisamp = (*thishit).energy();
0178           // find intercalib constant for this xtal
0179           auto icalit2 = icalMap.find((*thishit).id());
0180           EcalIntercalibConstant icalconst2 = 1.;
0181           if (icalit2 != icalMap.end()) {
0182             icalconst2 = (*icalit2);
0183             //     LogDebug("EcalRecHitDebug") << "Found intercalib for xtal " << EBDetId(it->id()).ic() << " " << icalconst ;
0184           } else {
0185             //edm::LogError("EcalRecHitError") << "No intercalib const found for xtal " << EBDetId(ebDet) << "! something wrong with EcalIntercalibConstants in your DB? " ;
0186           }
0187           float lasercalib2 = pLaser->getLaserCorrection(EBDetId((*thishit).id()), iEvent.time());
0188           thisamp /= (icalconst2 * lasercalib2 * adcconst);  ///LASER and CALIB constants from the DB
0189           if (thisamp > secondMin)
0190             secondMin = thisamp;
0191         }
0192       }
0193
0194       if (secondMin > minAmp2_) {
0195         thereIsSignal = true;
0196         break;
0197       }
0198     }
0199   }
0200   //std::cout << " Ok is There one of THEM " << thereIsSignal << std::endl;
0201   return thereIsSignal;
0202 }
0203
0204 void EcalMIPRecHitFilter::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0205   edm::ParameterSetDescription desc;
0206
0207   desc.add<edm::InputTag>("EcalRecHitCollection", edm::InputTag("ecalRecHit", "EcalRecHitsEB"));
0208   desc.addUntracked<double>("AmpMinSeed", 0.045);
0209   desc.addUntracked<double>("AmpMin2", 0.045);
0210   desc.addUntracked<double>("SingleAmpMin", 0.108);
0211   desc.addUntracked<std::vector<int>>("maskedChannels", std::vector<int>{});
0212   desc.addUntracked<int>("side", 3);
0213
0214   descriptions.add("ecalMIPRecHitFilter", desc);
0215 }
0216
0217 // declare this class as a framework plugin
0218 #include "FWCore/Framework/interface/MakerMacros.h"
0219 DEFINE_FWK_MODULE(EcalMIPRecHitFilter);