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