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     //number of accepted events
0192     unsigned int naccept = 0;
0193
0194     for (unsigned int itry = 0; itry < nAttempts_; ++itry) {
0195       hadronizer_.setLHEEvent(std::make_unique<lhef::LHEEvent>(hadronizer_.getLHERunInfo(), *product));
0196
0197       // hadronizer_.generatePartons();
0198       if (!hadronizer_.hadronize())
0199         continue;
0200
0201       //  this is "fake" stuff
0202       // in principle, decays are done as part of full event generation,
0203       // except for particles that are marked as to be kept stable
0204       // but we currently keep in it the design, because we might want
0205       // to use such feature for other applications
0206       //
0207       if (!hadronizer_.decay())
0208         continue;
0209
0210       std::unique_ptr<HepMC::GenEvent> event(hadronizer_.getGenEvent());
0211       if (!event.get())
0212         continue;
0213
0214       // The external decay driver is being added to the system,
0215       // it should be called here
0216       //
0217       if (decayer_) {
0218         auto lheEvent = hadronizer_.getLHEEvent();
0219         auto t = decayer_->decay(event.get(), lheEvent.get());
0220         if (t != event.get()) {
0221           event.reset(t);
0222         }
0223         hadronizer_.setLHEEvent(std::move(lheEvent));
0224       }
0225
0226       if (!event.get())
0227         continue;
0228
0229       // check and perform if there're any unstable particles after
0230       // running external decay packges
0231       //
0232       hadronizer_.resetEvent(std::move(event));
0233       if (!hadronizer_.residualDecay())
0234         continue;
0235
0236       hadronizer_.finalizeEvent();
0237
0238       event = hadronizer_.getGenEvent();
0239       if (!event.get())
0240         continue;
0241
0242       event->set_event_number(ev.id().event());
0243
0244       std::unique_ptr<GenEventInfoProduct> genEventInfo(hadronizer_.getGenEventInfo());
0245       if (!genEventInfo.get()) {
0246         // create GenEventInfoProduct from HepMC event in case hadronizer didn't provide one
0247         genEventInfo = std::make_unique<GenEventInfoProduct>(event.get());
0248       }
0249
0250       //if HepMCFilter was specified, test event
0251       if (filter_ && !filter_->filter(event.get(), genEventInfo->weight()))
0252         continue;
0253
0254       ++naccept;
0255
0256       //keep the LAST accepted event (which is equivalent to choosing randomly from the accepted events)
0257       finalEvent = std::move(event);
0258       finalGenEventInfo = std::move(genEventInfo);
0259     }
0260
0261     if (!naccept)
0262       return false;
0263
0264     //adjust event weights if necessary (in case input event was attempted multiple times)
0265     if (nAttempts_ > 1) {
0266       double multihadweight = double(naccept) / double(nAttempts_);
0267
0268       //adjust weight for GenEventInfoProduct
0269       finalGenEventInfo->weights()[0] *= multihadweight;
0270
0271       //adjust weight for HepMC GenEvent (used e.g for RIVET)
0272       finalEvent->weights()[0] *= multihadweight;
0273     }
0274
0275     ev.put(std::move(finalGenEventInfo));
0276
0277     std::unique_ptr<HepMCProduct> bare_product(new HepMCProduct());
0278     bare_product->addHepMCData(finalEvent.release());
0279     ev.put(std::move(bare_product), "unsmeared");
0280
0281     return true;
0282   }
0283
0284   template <class HAD, class DEC>
0285   void HadronizerFilter<HAD, DEC>::beginRun(Run const& run, EventSetup const& es) {
0286     // this is run-specific
0287
0288     // get LHE stuff and pass to hadronizer!
0289
0290     if (counterRunInfoProducts_ > 1)
0291       throw edm::Exception(errors::EventCorruption) << "More than one LHERunInfoProduct present";
0292
0293     if (counterRunInfoProducts_ == 0)
0294       throw edm::Exception(errors::EventCorruption) << "No LHERunInfoProduct present";
0295
0296     edm::Handle<LHERunInfoProduct> lheRunInfoProduct;
0297     run.getByLabel(runInfoProductTag_, lheRunInfoProduct);
0298     //TODO: fix so that this actually works with getByToken commented below...
0299     //run.getByToken(runInfoProductToken_, lheRunInfoProduct);
0300
0301     hadronizer_.setLHERunInfo(std::make_unique<lhef::LHERunInfo>(*lheRunInfoProduct));
0302     lhef::LHERunInfo* lheRunInfo = hadronizer_.getLHERunInfo().get();
0303     lheRunInfo->initLumi();
0304   }
0305
0306   template <class HAD, class DEC>
0307   void HadronizerFilter<HAD, DEC>::endRun(Run const& r, EventSetup const&) {}
0308
0309   template <class HAD, class DEC>
0310   void HadronizerFilter<HAD, DEC>::endRunProduce(Run& r, EventSetup const&) {
0311     // Retrieve the LHE run info summary and transfer determined
0312     // cross-section into the generator run info
0313
0314     const lhef::LHERunInfo* lheRunInfo = hadronizer_.getLHERunInfo().get();
0315     lhef::LHERunInfo::XSec xsec = lheRunInfo->xsec();
0316
0317     GenRunInfoProduct& genRunInfo = hadronizer_.getGenRunInfo();
0318     genRunInfo.setInternalXSec(GenRunInfoProduct::XSec(xsec.value(), xsec.error()));
0319
0320     // If relevant, record the integrated luminosity for this run
0321     // here.  To do so, we would need a standard function to invoke on
0322     // the contained hadronizer that would report the integrated
0323     // luminosity.
0324
0325     hadronizer_.statistics();
0326     if (decayer_)
0327       decayer_->statistics();
0328     if (filter_)
0329       filter_->statistics();
0330     lheRunInfo->statistics();
0331
0332     std::unique_ptr<GenRunInfoProduct> griproduct(new GenRunInfoProduct(genRunInfo));
0333     r.put(std::move(griproduct));
0334   }
0335
0336   template <class HAD, class DEC>
0337   void HadronizerFilter<HAD, DEC>::beginLuminosityBlock(LuminosityBlock const& lumi, EventSetup const& es) {}
0338
0339   template <class HAD, class DEC>
0340   void HadronizerFilter<HAD, DEC>::beginLuminosityBlockProduce(LuminosityBlock& lumi, EventSetup const& es) {
0341     lhef::LHERunInfo* lheRunInfo = hadronizer_.getLHERunInfo().get();
0342     lheRunInfo->initLumi();
0343
0344     RandomEngineSentry<HAD> randomEngineSentry(&hadronizer_, lumi.index());
0345     RandomEngineSentry<DEC> randomEngineSentryDecay(decayer_, lumi.index());
0346
0347     hadronizer_.randomizeIndex(lumi, randomEngineSentry.randomEngine());
0348
0349     if (!hadronizer_.readSettings(1))
0350       throw edm::Exception(errors::Configuration)
0351           << "Failed to read settings for the hadronizer " << hadronizer_.classname() << " \n";
0352
0353     if (decayer_) {
0354       decayer_->init(es);
0355       if (!hadronizer_.declareStableParticles(decayer_->operatesOnParticles()))
0356         throw edm::Exception(errors::Configuration) << "Failed to declare stable particles in hadronizer "
0357                                                     << hadronizer_.classname() << " for internal parton generation\n";
0358       if (!hadronizer_.declareSpecialSettings(decayer_->specialSettings()))
0359         throw edm::Exception(errors::Configuration)
0360             << "Failed to declare special settings in hadronizer " << hadronizer_.classname() << "\n";
0361     }
0362
0363     if (filter_) {
0364       filter_->resetStatistics();
0365     }
0366
0367     if (!hadronizer_.initializeForExternalPartons())
0368       throw edm::Exception(errors::Configuration)
0369           << "Failed to initialize hadronizer " << hadronizer_.classname() << " for external parton generation\n";
0370
0371     std::unique_ptr<GenLumiInfoHeader> genLumiInfoHeader(hadronizer_.getGenLumiInfoHeader());
0372     lumi.put(std::move(genLumiInfoHeader));
0373   }
0374
0375   template <class HAD, class DEC>
0376   void HadronizerFilter<HAD, DEC>::endLuminosityBlock(LuminosityBlock const&, EventSetup const&) {}
0377
0378   template <class HAD, class DEC>
0379   void HadronizerFilter<HAD, DEC>::endLuminosityBlockProduce(LuminosityBlock& lumi, EventSetup const&) {
0380     const lhef::LHERunInfo* lheRunInfo = hadronizer_.getLHERunInfo().get();
0381
0382     std::vector<lhef::LHERunInfo::Process> LHELumiProcess = lheRunInfo->getLumiProcesses();
0383     std::vector<GenLumiInfoProduct::ProcessInfo> GenLumiProcess;
0384     for (unsigned int i = 0; i < LHELumiProcess.size(); i++) {
0385       lhef::LHERunInfo::Process thisProcess = LHELumiProcess[i];
0386
0387       GenLumiInfoProduct::ProcessInfo temp;
0388       temp.setProcess(thisProcess.process());
0389       temp.setLheXSec(thisProcess.getLHEXSec().value(), thisProcess.getLHEXSec().error());
0390       temp.setNPassPos(thisProcess.nPassPos());
0391       temp.setNPassNeg(thisProcess.nPassNeg());
0392       temp.setNTotalPos(thisProcess.nTotalPos());
0393       temp.setNTotalNeg(thisProcess.nTotalNeg());
0394       temp.setTried(thisProcess.tried().n(), thisProcess.tried().sum(), thisProcess.tried().sum2());
0395       temp.setSelected(thisProcess.selected().n(), thisProcess.selected().sum(), thisProcess.selected().sum2());
0396       temp.setKilled(thisProcess.killed().n(), thisProcess.killed().sum(), thisProcess.killed().sum2());
0397       temp.setAccepted(thisProcess.accepted().n(), thisProcess.accepted().sum(), thisProcess.accepted().sum2());
0398       temp.setAcceptedBr(thisProcess.acceptedBr().n(), thisProcess.acceptedBr().sum(), thisProcess.acceptedBr().sum2());
0399       GenLumiProcess.push_back(temp);
0400     }
0401     std::unique_ptr<GenLumiInfoProduct> genLumiInfo(new GenLumiInfoProduct());
0402     genLumiInfo->setHEPIDWTUP(lheRunInfo->getHEPRUP()->IDWTUP);
0403     genLumiInfo->setProcessInfo(GenLumiProcess);
0404
0405     lumi.put(std::move(genLumiInfo));
0406
0407     // produce GenFilterInfo if HepMCFilter is called
0408     if (filter_) {
0409       std::unique_ptr<GenFilterInfo> thisProduct(new GenFilterInfo(filter_->numEventsPassPos(),
0410                                                                    filter_->numEventsPassNeg(),
0411                                                                    filter_->numEventsTotalPos(),
0412                                                                    filter_->numEventsTotalNeg(),
0413                                                                    filter_->sumpass_w(),
0414                                                                    filter_->sumpass_w2(),
0415                                                                    filter_->sumtotal_w(),
0416                                                                    filter_->sumtotal_w2()));
0417       lumi.put(std::move(thisProduct));
0418     }
0419   }
0420
0421 }  // namespace edm
0422
0423 #endif  // gen_HadronizerFilter_h