Core/plugins/GenFilterEfficiencyProducer.cc

0001 // F. Cossutti
0002 //
0003
0004 // producer of a summary information product on filter efficiency for a user specified path
0005 // meant for the generator filter efficiency calculation
0006
0007 // system include files
0008 #include <memory>
0009 #include <string>
0010 #include <atomic>
0011
0012 // user include files
0013 #include "FWCore/Framework/interface/Frameworkfwd.h"
0014 #include "FWCore/Framework/interface/global/EDProducer.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/Framework/interface/EventSetup.h"
0017 #include "FWCore/Framework/interface/Event.h"
0018 #include "FWCore/Framework/interface/LuminosityBlock.h"
0019 #include "FWCore/Framework/interface/MakerMacros.h"
0020 #include "FWCore/Utilities/interface/BranchType.h"
0021 #include "FWCore/Utilities/interface/EDGetToken.h"
0022
0023 #include "FWCore/ServiceRegistry/interface/Service.h"
0024 #include "FWCore/Framework/interface/TriggerNamesService.h"
0025
0026 #include "DataFormats/Common/interface/TriggerResults.h"
0027
0028 #include "SimDataFormats/GeneratorProducts/interface/GenFilterInfo.h"
0029 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0030 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0031
0032 namespace genFilterEff {
0033   struct Sums {
0034     mutable std::atomic<unsigned int> numEventsPassPos_ = {0};
0035     mutable std::atomic<unsigned int> numEventsPassNeg_ = {0};
0036     mutable std::atomic<unsigned int> numEventsTotalPos_ = {0};
0037     mutable std::atomic<unsigned int> numEventsTotalNeg_ = {0};
0038     mutable std::atomic<double> sumpass_w_ = {0};
0039     mutable std::atomic<double> sumpass_w2_ = {0};
0040     mutable std::atomic<double> sumtotal_w_ = {0};
0041     mutable std::atomic<double> sumtotal_w2_ = {0};
0042   };
0043 }  // namespace genFilterEff
0044
0045 namespace {
0046   void atomic_sum_double(std::atomic<double>& oValue, double element) {
0047     double v = oValue.load();
0048     double sum = v + element;
0049     while (not oValue.compare_exchange_strong(v, sum)) {
0050       //some other thread updated oValue
0051       sum = v + element;
0052     }
0053   }
0054 }  // namespace
0055 using namespace genFilterEff;
0056
0057 class GenFilterEfficiencyProducer
0058     : public edm::global::EDProducer<edm::EndLuminosityBlockProducer, edm::LuminosityBlockCache<Sums>> {
0059 public:
0060   explicit GenFilterEfficiencyProducer(const edm::ParameterSet&);
0061   ~GenFilterEfficiencyProducer() override;
0062
0063 private:
0064   void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0065
0066   void globalEndLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) const override;
0067   void globalEndLuminosityBlockProduce(edm::LuminosityBlock&, const edm::EventSetup&) const override;
0068
0069   std::shared_ptr<Sums> globalBeginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) const override;
0070   // ----------member data ---------------------------
0071
0072   edm::EDGetTokenT<edm::TriggerResults> triggerResultsToken_;
0073   edm::EDGetTokenT<GenEventInfoProduct> genEventInfoToken_;
0074
0075   std::string filterPath;
0076
0077   edm::service::TriggerNamesService* tns_;
0078
0079   std::string thisProcess;
0080   unsigned int pathIndex;
0081 };
0082
0083 GenFilterEfficiencyProducer::GenFilterEfficiencyProducer(const edm::ParameterSet& iConfig)
0084     : filterPath(iConfig.getParameter<std::string>("filterPath")), tns_(), thisProcess(), pathIndex(100000) {
0085   //now do what ever initialization is needed
0086   if (edm::Service<edm::service::TriggerNamesService>().isAvailable()) {
0087     // get tns pointer
0088     tns_ = edm::Service<edm::service::TriggerNamesService>().operator->();
0089     if (tns_ != nullptr) {
0090       thisProcess = tns_->getProcessName();
0091       std::vector<std::string> theNames = tns_->getTrigPaths();
0092       for (unsigned int i = 0; i < theNames.size(); i++) {
0093         if (theNames[i] == filterPath) {
0094           pathIndex = i;
0095           continue;
0096         }
0097       }
0098     } else
0099       edm::LogError("ServiceNotAvailable") << "TriggerNamesServive not available, no filter information stored";
0100   }
0101
0102   triggerResultsToken_ = consumes<edm::TriggerResults>(edm::InputTag("TriggerResults", "", thisProcess));
0103   genEventInfoToken_ = consumes<GenEventInfoProduct>(edm::InputTag("generator", ""));
0104   produces<GenFilterInfo, edm::Transition::EndLuminosityBlock>();
0105 }
0106
0107 GenFilterEfficiencyProducer::~GenFilterEfficiencyProducer() {}
0108
0109 //
0110 // member functions
0111 //
0112
0113 // ------------ method called to for each event  ------------
0114 void GenFilterEfficiencyProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0115   edm::Handle<edm::TriggerResults> trigR;
0116   iEvent.getByToken(triggerResultsToken_, trigR);
0117   edm::Handle<GenEventInfoProduct> genEventScale;
0118   iEvent.getByToken(genEventInfoToken_, genEventScale);
0119   if (!genEventScale.isValid())
0120     return;
0121   double weight = genEventScale->weight();
0122
0123   auto sums = luminosityBlockCache(iEvent.getLuminosityBlock().index());
0124
0125   unsigned int nSize = (*trigR).size();
0126   // std::cout << "Number of paths in TriggerResults = " << nSize  << std::endl;
0127   if (nSize >= pathIndex) {
0128     if (!trigR->wasrun(pathIndex))
0129       return;
0130     if (trigR->accept(pathIndex)) {
0131       atomic_sum_double(sums->sumpass_w_, weight);
0132       atomic_sum_double(sums->sumpass_w2_, weight * weight);
0133
0134       atomic_sum_double(sums->sumtotal_w_, weight);
0135       atomic_sum_double(sums->sumtotal_w2_, weight * weight);
0136
0137       if (weight > 0) {
0138         sums->numEventsPassPos_++;
0139         sums->numEventsTotalPos_++;
0140       } else {
0141         sums->numEventsPassNeg_++;
0142         sums->numEventsTotalNeg_++;
0143       }
0144
0145     } else  // if fail the filter
0146     {
0147       atomic_sum_double(sums->sumtotal_w_, weight);
0148       atomic_sum_double(sums->sumtotal_w2_, weight * weight);
0149
0150       if (weight > 0)
0151         sums->numEventsTotalPos_++;
0152       else
0153         sums->numEventsTotalNeg_++;
0154     }
0155     //    std::cout << "Total events = " << numEventsTotal << " passed = " << numEventsPassed << std::endl;
0156   }
0157 }
0158
0159 std::shared_ptr<Sums> GenFilterEfficiencyProducer::globalBeginLuminosityBlock(edm::LuminosityBlock const&,
0160                                                                               edm::EventSetup const&) const {
0161   return std::make_shared<Sums>();
0162 }
0163
0164 void GenFilterEfficiencyProducer::globalEndLuminosityBlock(edm::LuminosityBlock const& iLumi,
0165                                                            const edm::EventSetup&) const {}
0166
0167 void GenFilterEfficiencyProducer::globalEndLuminosityBlockProduce(edm::LuminosityBlock& iLumi,
0168                                                                   const edm::EventSetup&) const {
0169   auto sums = luminosityBlockCache(iLumi.index());
0170   std::unique_ptr<GenFilterInfo> thisProduct(new GenFilterInfo(sums->numEventsPassPos_,
0171                                                                sums->numEventsPassNeg_,
0172                                                                sums->numEventsTotalPos_,
0173                                                                sums->numEventsTotalNeg_,
0174                                                                sums->sumpass_w_,
0175                                                                sums->sumpass_w2_,
0176                                                                sums->sumtotal_w_,
0177                                                                sums->sumtotal_w2_));
0178   iLumi.put(std::move(thisProduct));
0179 }
0180
0181 DEFINE_FWK_MODULE(GenFilterEfficiencyProducer);