Core/interface/HadronizerFilter.h

0001 // -*- C++ -*-
0002 //
0003 //
0004
0005 // class template HadronizerFilter<HAD> provides an EDFilter which uses
0006 // the hadronizer type HAD to read in external partons and hadronize them,
0007 // and decay the resulting particles, in the CMS framework.
0008
0009 #ifndef gen_HadronizerFilter_h
0010 #define gen_HadronizerFilter_h
0011
0012 #include <memory>
0013 #include <string>
0014 #include <vector>
0015
0016 #include "FWCore/Concurrency/interface/SharedResourceNames.h"
0017 #include "FWCore/Framework/interface/one/EDFilter.h"
0018 #include "FWCore/Framework/interface/Event.h"
0019 #include "FWCore/Framework/interface/EventSetup.h"
0020 #include "FWCore/Framework/interface/FileBlock.h"
0021 #include "FWCore/Framework/interface/LuminosityBlock.h"
0022 #include "FWCore/Framework/interface/MakerMacros.h"
0023 #include "FWCore/Framework/interface/Run.h"
0024 #include "FWCore/Framework/interface/ConsumesCollector.h"
0025 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0026 #include "FWCore/ServiceRegistry/interface/RandomEngineSentry.h"
0027 #include "FWCore/Utilities/interface/BranchType.h"
0028 #include "FWCore/Utilities/interface/EDMException.h"
0029 #include "FWCore/Utilities/interface/EDGetToken.h"
0030 #include "FWCore/Utilities/interface/TypeID.h"
0031 #include "DataFormats/Provenance/interface/BranchDescription.h"
0032 #include "CLHEP/Random/RandomEngine.h"
0033
0034 // #include "GeneratorInterface/ExternalDecays/interface/ExternalDecayDriver.h"
0035
0036 #include "GeneratorInterface/Core/interface/HepMCFilterDriver.h"
0037
0038 // LHE Run
0039 #include "SimDataFormats/GeneratorProducts/interface/LHERunInfoProduct.h"
0040 #include "GeneratorInterface/LHEInterface/interface/LHERunInfo.h"
0041
0042 // LHE Event
0043 #include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h"
0044 #include "GeneratorInterface/LHEInterface/interface/LHEEvent.h"
0045
0046 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0047 #include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
0048 #include "SimDataFormats/GeneratorProducts/interface/GenLumiInfoHeader.h"
0049 #include "SimDataFormats/GeneratorProducts/interface/GenLumiInfoProduct.h"
0050 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0051
0052 namespace edm {
0053   template <class HAD, class DEC>
0054   class HadronizerFilter : public one::EDFilter<EndRunProducer,
0055                                                 BeginLuminosityBlockProducer,
0056                                                 EndLuminosityBlockProducer,
0057                                                 one::WatchRuns,
0058                                                 one::WatchLuminosityBlocks,
0059                                                 one::SharedResources> {
0060   public:
0061     typedef HAD Hadronizer;
0062     typedef DEC Decayer;
0063
0064     // The given ParameterSet will be passed to the contained
0065     // Hadronizer object.
0066     explicit HadronizerFilter(ParameterSet const& ps);
0067
0068     ~HadronizerFilter() override;
0069
0070     bool filter(Event& e, EventSetup const& es) override;
0071     void beginRun(Run const&, EventSetup const&) override;
0072     void endRun(Run const&, EventSetup const&) override;
0073     void endRunProduce(Run&, EventSetup const&) override;
0074     void beginLuminosityBlock(LuminosityBlock const&, EventSetup const&) override;
0075     void beginLuminosityBlockProduce(LuminosityBlock&, EventSetup const&) override;
0076     void endLuminosityBlock(LuminosityBlock const&, EventSetup const&) override;
0077     void endLuminosityBlockProduce(LuminosityBlock&, EventSetup const&) override;
0078
0079   private:
0080     Hadronizer hadronizer_;
0081     // gen::ExternalDecayDriver* decayer_;
0082     Decayer* decayer_;
0083     HepMCFilterDriver* filter_;
0084     InputTag runInfoProductTag_;
0085     EDGetTokenT<LHERunInfoProduct> runInfoProductToken_;
0086     EDGetTokenT<LHEEventProduct> eventProductToken_;
0087     unsigned int counterRunInfoProducts_;
0088     unsigned int nAttempts_;
0089   };
0090
0091   //------------------------------------------------------------------------
0092   //
0093   // Implementation
0094
0095   template <class HAD, class DEC>
0096   HadronizerFilter<HAD, DEC>::HadronizerFilter(ParameterSet const& ps)
0097       : EDFilter(),
0098         hadronizer_(ps),
0099         decayer_(nullptr),
0100         filter_(nullptr),
0101         runInfoProductTag_(),
0102         runInfoProductToken_(),
0103         eventProductToken_(),
0104         counterRunInfoProducts_(0),
0105         nAttempts_(1) {
0106     callWhenNewProductsRegistered([this](BranchDescription const& iBD) {
0107       //this is called each time a module registers that it will produce a LHERunInfoProduct
0108       if (iBD.unwrappedTypeID() == edm::TypeID(typeid(LHERunInfoProduct)) && iBD.branchType() == InRun) {
0109         ++(this->counterRunInfoProducts_);
0110         this->eventProductToken_ = consumes<LHEEventProduct>(
0111             InputTag((iBD.moduleLabel() == "externalLHEProducer") ? "externalLHEProducer" : "source"));
0112         this->runInfoProductTag_ = InputTag(iBD.moduleLabel(), iBD.productInstanceName(), iBD.processName());
0113         this->runInfoProductToken_ = consumes<LHERunInfoProduct, InRun>(
0114             InputTag(iBD.moduleLabel(), iBD.productInstanceName(), iBD.processName()));
0115       }
0116     });
0117
0118     // TODO:
0119     // Put the list of types produced by the filters here.
0120     // The current design calls for:
0121     //   * LHEGeneratorInfo
0122     //   * LHEEvent
0123     //   * HepMCProduct
0124     // But I can not find the LHEGeneratorInfo class; it might need to
0125     // be invented.
0126
0127     std::vector<std::string> const& sharedResources = hadronizer_.sharedResources();
0128     for (auto const& resource : sharedResources) {
0129       usesResource(resource);
0130     }
0131
0132     if (ps.exists("ExternalDecays")) {
0133       //decayer_ = new gen::ExternalDecayDriver(ps.getParameter<ParameterSet>("ExternalDecays"));
0134       ParameterSet ps1 = ps.getParameter<ParameterSet>("ExternalDecays");
0135       decayer_ = new Decayer(ps1, consumesCollector());
0136
0137       std::vector<std::string> const& sharedResourcesDec = decayer_->sharedResources();
0138       for (auto const& resource : sharedResourcesDec) {
0139         usesResource(resource);
0140       }
0141     }
0142
0143     if (ps.exists("HepMCFilter")) {
0144       ParameterSet psfilter = ps.getParameter<ParameterSet>("HepMCFilter");
0145       filter_ = new HepMCFilterDriver(psfilter);
0146     }
0147
0148     //initialize setting for multiple hadronization attempts
0149     if (ps.exists("nAttempts")) {
0150       nAttempts_ = ps.getParameter<unsigned int>("nAttempts");
0151     }
0152
0153     // This handles the case where there are no shared resources, because you
0154     // have to declare something when the SharedResources template parameter was used.
0155     if (sharedResources.empty() && (!decayer_ || decayer_->sharedResources().empty())) {
0156       usesResource(edm::uniqueSharedResourceName());
0157     }
0158
0159     produces<edm::HepMCProduct>("unsmeared");
0160     produces<GenEventInfoProduct>();
0161     produces<GenLumiInfoHeader, edm::Transition::BeginLuminosityBlock>();
0162     produces<GenLumiInfoProduct, edm::Transition::EndLuminosityBlock>();
0163     produces<GenRunInfoProduct, edm::Transition::EndRun>();
0164     if (filter_)
0165       produces<GenFilterInfo, edm::Transition::EndLuminosityBlock>();
0166   }
0167
0168   template <class HAD, class DEC>
0169   HadronizerFilter<HAD, DEC>::~HadronizerFilter() {
0170     if (decayer_)
0171       delete decayer_;
0172     if (filter_)
0173       delete filter_;
0174   }
0175
0176   template <class HAD, class DEC>
0177   bool HadronizerFilter<HAD, DEC>::filter(Event& ev, EventSetup const& /* es */) {
0178     RandomEngineSentry<HAD> randomEngineSentry(&hadronizer_, ev.streamID());
0179     RandomEngineSentry<DEC> randomEngineSentryDecay(decayer_, ev.streamID());
0180
0181     hadronizer_.setEDMEvent(ev);
0182
0183     // get LHE stuff and pass to hadronizer!
0184     //
0185     edm::Handle<LHEEventProduct> product;
0186     ev.getByToken(eventProductToken_, product);
0187
0188     std::unique_ptr<HepMC::GenEvent> finalEvent;
0189     std::unique_ptr<GenEventInfoProduct> finalGenEventInfo;
0190
0191     //sum of weights for events passing hadronization
0192     double waccept = 0;
0193     //number of accepted events
0194     unsigned int naccept = 0;
0195
0196     for (unsigned int itry = 0; itry < nAttempts_; ++itry) {
0197       hadronizer_.setLHEEvent(std::make_unique<lhef::LHEEvent>(hadronizer_.getLHERunInfo(), *product));
0198
0199       // hadronizer_.generatePartons();
0200       if (!hadronizer_.hadronize())
0201         continue;
0202
0203       //  this is "fake" stuff
0204       // in principle, decays are done as part of full event generation,
0205       // except for particles that are marked as to be kept stable
0206       // but we currently keep in it the design, because we might want
0207       // to use such feature for other applications
0208       //
0209       if (!hadronizer_.decay())
0210         continue;
0211
0212       std::unique_ptr<HepMC::GenEvent> event(hadronizer_.getGenEvent());
0213       if (!event.get())
0214         continue;
0215
0216       // The external decay driver is being added to the system,
0217       // it should be called here
0218       //
0219       if (decayer_) {
0220         auto lheEvent = hadronizer_.getLHEEvent();
0221         auto t = decayer_->decay(event.get(), lheEvent.get());
0222         if (t != event.get()) {
0223           event.reset(t);
0224         }
0225         hadronizer_.setLHEEvent(std::move(lheEvent));
0226       }
0227
0228       if (!event.get())
0229         continue;
0230
0231       // check and perform if there're any unstable particles after
0232       // running external decay packges
0233       //
0234       hadronizer_.resetEvent(std::move(event));
0235       if (!hadronizer_.residualDecay())
0236         continue;
0237
0238       hadronizer_.finalizeEvent();
0239
0240       event = hadronizer_.getGenEvent();
0241       if (!event.get())
0242         continue;
0243
0244       event->set_event_number(ev.id().event());
0245
0246       std::unique_ptr<GenEventInfoProduct> genEventInfo(hadronizer_.getGenEventInfo());
0247       if (!genEventInfo.get()) {
0248         // create GenEventInfoProduct from HepMC event in case hadronizer didn't provide one
0249         genEventInfo = std::make_unique<GenEventInfoProduct>(event.get());
0250       }
0251
0252       //if HepMCFilter was specified, test event
0253       if (filter_ && !filter_->filter(event.get(), genEventInfo->weight()))
0254         continue;
0255
0256       waccept += genEventInfo->weight();
0257       ++naccept;
0258
0259       //keep the LAST accepted event (which is equivalent to choosing randomly from the accepted events)
0260       finalEvent = std::move(event);
0261       finalGenEventInfo = std::move(genEventInfo);
0262     }
0263
0264     if (!naccept)
0265       return false;
0266
0267     //adjust event weights if necessary (in case input event was attempted multiple times)
0268     if (nAttempts_ > 1) {
0269       double multihadweight = double(naccept) / double(nAttempts_);
0270
0271       //adjust weight for GenEventInfoProduct
0272       finalGenEventInfo->weights()[0] *= multihadweight;
0273
0274       //adjust weight for HepMC GenEvent (used e.g for RIVET)
0275       finalEvent->weights()[0] *= multihadweight;
0276     }
0277
0278     ev.put(std::move(finalGenEventInfo));
0279
0280     std::unique_ptr<HepMCProduct> bare_product(new HepMCProduct());
0281     bare_product->addHepMCData(finalEvent.release());
0282     ev.put(std::move(bare_product), "unsmeared");
0283
0284     return true;
0285   }
0286
0287   template <class HAD, class DEC>
0288   void HadronizerFilter<HAD, DEC>::beginRun(Run const& run, EventSetup const& es) {
0289     // this is run-specific
0290
0291     // get LHE stuff and pass to hadronizer!
0292
0293     if (counterRunInfoProducts_ > 1)
0294       throw edm::Exception(errors::EventCorruption) << "More than one LHERunInfoProduct present";
0295
0296     if (counterRunInfoProducts_ == 0)
0297       throw edm::Exception(errors::EventCorruption) << "No LHERunInfoProduct present";
0298
0299     edm::Handle<LHERunInfoProduct> lheRunInfoProduct;
0300     run.getByLabel(runInfoProductTag_, lheRunInfoProduct);
0301     //TODO: fix so that this actually works with getByToken commented below...
0302     //run.getByToken(runInfoProductToken_, lheRunInfoProduct);
0303
0304     hadronizer_.setLHERunInfo(std::make_unique<lhef::LHERunInfo>(*lheRunInfoProduct));
0305     lhef::LHERunInfo* lheRunInfo = hadronizer_.getLHERunInfo().get();
0306     lheRunInfo->initLumi();
0307   }
0308
0309   template <class HAD, class DEC>
0310   void HadronizerFilter<HAD, DEC>::endRun(Run const& r, EventSetup const&) {}
0311
0312   template <class HAD, class DEC>
0313   void HadronizerFilter<HAD, DEC>::endRunProduce(Run& r, EventSetup const&) {
0314     // Retrieve the LHE run info summary and transfer determined
0315     // cross-section into the generator run info
0316
0317     const lhef::LHERunInfo* lheRunInfo = hadronizer_.getLHERunInfo().get();
0318     lhef::LHERunInfo::XSec xsec = lheRunInfo->xsec();
0319
0320     GenRunInfoProduct& genRunInfo = hadronizer_.getGenRunInfo();
0321     genRunInfo.setInternalXSec(GenRunInfoProduct::XSec(xsec.value(), xsec.error()));
0322
0323     // If relevant, record the integrated luminosity for this run
0324     // here.  To do so, we would need a standard function to invoke on
0325     // the contained hadronizer that would report the integrated
0326     // luminosity.
0327
0328     hadronizer_.statistics();
0329     if (decayer_)
0330       decayer_->statistics();
0331     if (filter_)
0332       filter_->statistics();
0333     lheRunInfo->statistics();
0334
0335     std::unique_ptr<GenRunInfoProduct> griproduct(new GenRunInfoProduct(genRunInfo));
0336     r.put(std::move(griproduct));
0337   }
0338
0339   template <class HAD, class DEC>
0340   void HadronizerFilter<HAD, DEC>::beginLuminosityBlock(LuminosityBlock const& lumi, EventSetup const& es) {}
0341
0342   template <class HAD, class DEC>
0343   void HadronizerFilter<HAD, DEC>::beginLuminosityBlockProduce(LuminosityBlock& lumi, EventSetup const& es) {
0344     lhef::LHERunInfo* lheRunInfo = hadronizer_.getLHERunInfo().get();
0345     lheRunInfo->initLumi();
0346
0347     RandomEngineSentry<HAD> randomEngineSentry(&hadronizer_, lumi.index());
0348     RandomEngineSentry<DEC> randomEngineSentryDecay(decayer_, lumi.index());
0349
0350     hadronizer_.randomizeIndex(lumi, randomEngineSentry.randomEngine());
0351
0352     if (!hadronizer_.readSettings(1))
0353       throw edm::Exception(errors::Configuration)
0354           << "Failed to read settings for the hadronizer " << hadronizer_.classname() << " \n";
0355
0356     if (decayer_) {
0357       decayer_->init(es);
0358       if (!hadronizer_.declareStableParticles(decayer_->operatesOnParticles()))
0359         throw edm::Exception(errors::Configuration) << "Failed to declare stable particles in hadronizer "
0360                                                     << hadronizer_.classname() << " for internal parton generation\n";
0361       if (!hadronizer_.declareSpecialSettings(decayer_->specialSettings()))
0362         throw edm::Exception(errors::Configuration)
0363             << "Failed to declare special settings in hadronizer " << hadronizer_.classname() << "\n";
0364     }
0365
0366     if (filter_) {
0367       filter_->resetStatistics();
0368     }
0369
0370     if (!hadronizer_.initializeForExternalPartons())
0371       throw edm::Exception(errors::Configuration)
0372           << "Failed to initialize hadronizer " << hadronizer_.classname() << " for external parton generation\n";
0373
0374     std::unique_ptr<GenLumiInfoHeader> genLumiInfoHeader(hadronizer_.getGenLumiInfoHeader());
0375     lumi.put(std::move(genLumiInfoHeader));
0376   }
0377
0378   template <class HAD, class DEC>
0379   void HadronizerFilter<HAD, DEC>::endLuminosityBlock(LuminosityBlock const&, EventSetup const&) {}
0380
0381   template <class HAD, class DEC>
0382   void HadronizerFilter<HAD, DEC>::endLuminosityBlockProduce(LuminosityBlock& lumi, EventSetup const&) {
0383     const lhef::LHERunInfo* lheRunInfo = hadronizer_.getLHERunInfo().get();
0384
0385     std::vector<lhef::LHERunInfo::Process> LHELumiProcess = lheRunInfo->getLumiProcesses();
0386     std::vector<GenLumiInfoProduct::ProcessInfo> GenLumiProcess;
0387     for (unsigned int i = 0; i < LHELumiProcess.size(); i++) {
0388       lhef::LHERunInfo::Process thisProcess = LHELumiProcess[i];
0389
0390       GenLumiInfoProduct::ProcessInfo temp;
0391       temp.setProcess(thisProcess.process());
0392       temp.setLheXSec(thisProcess.getLHEXSec().value(), thisProcess.getLHEXSec().error());
0393       temp.setNPassPos(thisProcess.nPassPos());
0394       temp.setNPassNeg(thisProcess.nPassNeg());
0395       temp.setNTotalPos(thisProcess.nTotalPos());
0396       temp.setNTotalNeg(thisProcess.nTotalNeg());
0397       temp.setTried(thisProcess.tried().n(), thisProcess.tried().sum(), thisProcess.tried().sum2());
0398       temp.setSelected(thisProcess.selected().n(), thisProcess.selected().sum(), thisProcess.selected().sum2());
0399       temp.setKilled(thisProcess.killed().n(), thisProcess.killed().sum(), thisProcess.killed().sum2());
0400       temp.setAccepted(thisProcess.accepted().n(), thisProcess.accepted().sum(), thisProcess.accepted().sum2());
0401       temp.setAcceptedBr(thisProcess.acceptedBr().n(), thisProcess.acceptedBr().sum(), thisProcess.acceptedBr().sum2());
0402       GenLumiProcess.push_back(temp);
0403     }
0404     std::unique_ptr<GenLumiInfoProduct> genLumiInfo(new GenLumiInfoProduct());
0405     genLumiInfo->setHEPIDWTUP(lheRunInfo->getHEPRUP()->IDWTUP);
0406     genLumiInfo->setProcessInfo(GenLumiProcess);
0407
0408     lumi.put(std::move(genLumiInfo));
0409
0410     // produce GenFilterInfo if HepMCFilter is called
0411     if (filter_) {
0412       std::unique_ptr<GenFilterInfo> thisProduct(new GenFilterInfo(filter_->numEventsPassPos(),
0413                                                                    filter_->numEventsPassNeg(),
0414                                                                    filter_->numEventsTotalPos(),
0415                                                                    filter_->numEventsTotalNeg(),
0416                                                                    filter_->sumpass_w(),
0417                                                                    filter_->sumpass_w2(),
0418                                                                    filter_->sumtotal_w(),
0419                                                                    filter_->sumtotal_w2()));
0420       lumi.put(std::move(thisProduct));
0421     }
0422   }
0423
0424 }  // namespace edm
0425
0426 #endif  // gen_HadronizerFilter_h