File indexing completed on 2025-03-23 15:57:48
0001
0002
0003
0004
0005
0006
0007
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
0035
0036 #include "GeneratorInterface/Core/interface/HepMCFilterDriver.h"
0037 #include "GeneratorInterface/Core/interface/HepMC3FilterDriver.h"
0038
0039
0040 #include "SimDataFormats/GeneratorProducts/interface/LHERunInfoProduct.h"
0041 #include "GeneratorInterface/LHEInterface/interface/LHERunInfo.h"
0042
0043
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
0068
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
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
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
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
0125
0126
0127
0128
0129
0130
0131
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
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
0160 if (ps.exists("nAttempts")) {
0161 nAttempts_ = ps.getParameter<unsigned int>("nAttempts");
0162 }
0163
0164
0165
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& ) {
0196 RandomEngineSentry<HAD> randomEngineSentry(&hadronizer_, ev.streamID());
0197 RandomEngineSentry<DEC> randomEngineSentryDecay(decayer_, ev.streamID());
0198
0199 hadronizer_.setEDMEvent(ev);
0200
0201
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
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
0218 if (!hadronizer_.hadronize())
0219 continue;
0220
0221
0222
0223
0224
0225
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
0238
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
0255
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) {
0272 event->set_event_number(ev.id().event());
0273 std::unique_ptr<GenEventInfoProduct> genEventInfo(hadronizer_.getGenEventInfo());
0274 if (!genEventInfo.get()) {
0275
0276 genEventInfo = std::make_unique<GenEventInfoProduct>(event.get());
0277 }
0278
0279
0280 if (filter_ && !filter_->filter(event.get(), genEventInfo->weight()))
0281 continue;
0282
0283 ++naccept;
0284
0285
0286 finalEvent = std::move(event);
0287 finalGenEventInfo = std::move(genEventInfo);
0288 } else if (ivhepmc == 3) {
0289 event3->set_event_number(ev.id().event());
0290 std::unique_ptr<GenEventInfoProduct3> genEventInfo3(hadronizer_.getGenEventInfo3());
0291 if (!genEventInfo3.get()) {
0292
0293 genEventInfo3 = std::make_unique<GenEventInfoProduct3>(event3.get());
0294 }
0295
0296
0297 if (filter3_ && !filter3_->filter(event3.get(), genEventInfo3->weight()))
0298 continue;
0299
0300 ++naccept;
0301
0302
0303 finalEvent3 = std::move(event3);
0304 finalGenEventInfo3 = std::move(genEventInfo3);
0305 }
0306 }
0307
0308 if (!naccept)
0309 return false;
0310
0311 if (ivhepmc == 2) {
0312
0313 if (nAttempts_ > 1) {
0314 double multihadweight = double(naccept) / double(nAttempts_);
0315
0316
0317 finalGenEventInfo->weights()[0] *= multihadweight;
0318
0319
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) {
0329
0330 if (nAttempts_ > 1) {
0331 double multihadweight = double(naccept) / double(nAttempts_);
0332
0333
0334 finalGenEventInfo3->weights()[0] *= multihadweight;
0335
0336
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
0353
0354
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
0365
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
0378
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
0387
0388
0389
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
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 }
0504
0505 #endif