HiGenCommon/plugins/HiCentralityBiasFilter.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    HiCentralityBiasFilter
0004 // Class:      HiCentralityBiasFilter
0005 //
0006 /**\class HiCentralityBiasFilter HiCentralityBiasFilter.cc yetkin/HiCentralityBiasFilter/src/HiCentralityBiasFilter.cc
0007
0008  Description: <one line class summary>
0009
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Yetkin Yilmaz
0015 //         Created:  Tue Aug 11 12:42:25 EDT 2009
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/one/EDFilter.h"
0025
0026 #include "FWCore/Framework/interface/Event.h"
0027 #include "FWCore/Framework/interface/MakerMacros.h"
0028
0029 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0030
0031 #include "FWCore/Utilities/interface/EDGetToken.h"
0032 #include "FWCore/Utilities/interface/InputTag.h"
0033 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0034
0035 #include "FWCore/ServiceRegistry/interface/Service.h"
0036
0037 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0038
0039 #include "CLHEP/Random/RandomEngine.h"
0040 #include "HepMC/GenEvent.h"
0041 #include "HepMC/HeavyIon.h"
0042 #include "TF1.h"
0043
0044 using namespace std;
0045
0046 //
0047 // class declaration
0048 //
0049
0050 class HiCentralityBiasFilter : public edm::one::EDFilter<> {
0051 public:
0052   explicit HiCentralityBiasFilter(const edm::ParameterSet&);
0053   ~HiCentralityBiasFilter() override = default;
0054
0055 private:
0056   void beginJob() override;
0057   bool filter(edm::Event&, const edm::EventSetup&) override;
0058   void endJob() override;
0059
0060   const edm::EDGetTokenT<edm::HepMCProduct> hepmcSrc_;
0061   const std::string func_;
0062   const std::vector<double> par_;
0063
0064   edm::Service<edm::RandomNumberGenerator> rng_;
0065
0066   TF1* fBias_;
0067 };
0068
0069 //
0070 // constants, enums and typedefs
0071 //
0072
0073 //
0074 // static data member definitions
0075 //
0076
0077 //
0078 // constructors and destructor
0079 //
0080 HiCentralityBiasFilter::HiCentralityBiasFilter(const edm::ParameterSet& iConfig)
0081     : hepmcSrc_(consumes<edm::HepMCProduct>(iConfig.getParameter<edm::InputTag>("generatorSmeared"))),
0082       func_(iConfig.getParameter<string>("function")),
0083       par_(iConfig.getParameter<vector<double> >("parameters")) {
0084   // Do what ever initialization is needed
0085 }
0086
0087 //
0088 // member functions
0089 //
0090
0091 // ------------ method called on each new Event  ------------
0092 bool HiCentralityBiasFilter::filter(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0093   using namespace edm;
0094
0095   CLHEP::HepRandomEngine& engine = rng_->getEngine(iEvent.streamID());
0096
0097   const edm::Handle<HepMCProduct>& mc = iEvent.getHandle(hepmcSrc_);
0098   const HepMC::GenEvent* evt = mc->GetEvent();
0099
0100   const HepMC::HeavyIon* hi = evt->heavy_ion();
0101   if (!hi)
0102     return false;
0103
0104   double b = hi->impact_parameter();
0105   double bound = fBias_->Eval(b);
0106   double rand = engine.flat();
0107   if (rand > bound)
0108     return false;
0109
0110   return true;
0111 }
0112
0113 // ------------ method called once each job just before starting event loop  ------------
0114 void HiCentralityBiasFilter::beginJob() {
0115   fBias_ = new TF1("fBias", func_.data(), 0, 20);
0116
0117   for (size_t ip = 0; ip < par_.size(); ++ip) {
0118     fBias_->SetParameter(ip, par_[ip]);
0119   }
0120
0121   double maxpoint = fBias_->GetMaximum(-0.1, 20);
0122   if (maxpoint < 0.9)
0123     throw cms::Exception("HeavyIonCentralityBias")
0124         << "Input bias function is not optimized. Peak value is " << maxpoint
0125         << " which is required to be close to 1. Please fix the parameters before production." << endl;
0126 }
0127
0128 // ------------ method called once each job just after ending the event loop  ------------
0129 void HiCentralityBiasFilter::endJob() {}
0130
0131 //define this as a plug-in
0132 DEFINE_FWK_MODULE(HiCentralityBiasFilter);