Core/interface/GeneratorFilter.h

0001 // -*- C++ -*-
0002 //
0003 //
0004
0005 // class template GeneratorFilter<HAD> provides an EDFilter which uses
0006 // the hadronizer type HAD to generate partons, hadronize them, and
0007 // decay the resulting particles, in the CMS framework.
0008
0009 #ifndef gen_GeneratorFilter_h
0010 #define gen_GeneratorFilter_h
0011
0012 #include <memory>
0013 #include <string>
0014 #include <vector>
0015
0016 #include "HepMC3/GenEvent.h"
0017
0018 #include "FWCore/Concurrency/interface/SharedResourceNames.h"
0019 #include "FWCore/Framework/interface/one/EDFilter.h"
0020 #include "FWCore/Framework/interface/Event.h"
0021 #include "FWCore/Framework/interface/EventSetup.h"
0022 #include "FWCore/Framework/interface/FileBlock.h"
0023 #include "FWCore/Framework/interface/LuminosityBlock.h"
0024 #include "FWCore/Framework/interface/MakerMacros.h"
0025 #include "FWCore/Framework/interface/Run.h"
0026 #include "FWCore/Framework/interface/ConsumesCollector.h"
0027 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0028 #include "FWCore/ServiceRegistry/interface/RandomEngineSentry.h"
0029 #include "FWCore/Utilities/interface/EDMException.h"
0030 #include "CLHEP/Random/RandomEngine.h"
0031
0032 // #include "GeneratorInterface/ExternalDecays/interface/ExternalDecayDriver.h"
0033
0034 //#include "GeneratorInterface/LHEInterface/interface/LHEEvent.h"
0035 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0036 #include "SimDataFormats/GeneratorProducts/interface/HepMC3Product.h"
0037 #include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
0038 #include "SimDataFormats/GeneratorProducts/interface/GenLumiInfoHeader.h"
0039 #include "SimDataFormats/GeneratorProducts/interface/GenLumiInfoProduct.h"
0040 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0041 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct3.h"
0042
0043 namespace edm {
0044   template <class HAD, class DEC>
0045   class GeneratorFilter : public one::EDFilter<EndRunProducer,
0046                                                BeginLuminosityBlockProducer,
0047                                                EndLuminosityBlockProducer,
0048                                                one::WatchLuminosityBlocks,
0049                                                one::SharedResources> {
0050   public:
0051     typedef HAD Hadronizer;
0052     typedef DEC Decayer;
0053
0054     // The given ParameterSet will be passed to the contained
0055     // Hadronizer object.
0056     explicit GeneratorFilter(ParameterSet const& ps);
0057
0058     ~GeneratorFilter() override;
0059
0060     bool filter(Event& e, EventSetup const& es) override;
0061     void endRunProduce(Run&, EventSetup const&) override;
0062     void beginLuminosityBlock(LuminosityBlock const&, EventSetup const&) override;
0063     void beginLuminosityBlockProduce(LuminosityBlock&, EventSetup const&) override;
0064     void endLuminosityBlock(LuminosityBlock const&, EventSetup const&) override;
0065     void endLuminosityBlockProduce(LuminosityBlock&, EventSetup const&) override;
0066     void preallocThreads(unsigned int iThreads) override;
0067
0068   private:
0069     // two-phase construction to allow specialization of the constructor
0070     void init(ParameterSet const& ps);
0071
0072     Hadronizer hadronizer_;
0073     //gen::ExternalDecayDriver* decayer_;
0074     Decayer* decayer_ = nullptr;
0075     unsigned int nEventsInLumiBlock_ = 0;
0076     unsigned int nThreads_{1};
0077     bool initialized_ = false;
0078     unsigned int ivhepmc = 2;
0079   };
0080
0081   //------------------------------------------------------------------------
0082   //
0083   // Implementation
0084
0085   template <class HAD, class DEC>
0086   GeneratorFilter<HAD, DEC>::GeneratorFilter(ParameterSet const& ps) : hadronizer_(ps) {
0087     init(ps);
0088   }
0089
0090   template <class HAD, class DEC>
0091   void GeneratorFilter<HAD, DEC>::init(ParameterSet const& ps) {
0092     // TODO:
0093     // Put the list of types produced by the filters here.
0094     // The current design calls for:
0095     //   * GenRunInfoProduct
0096     //   * HepMCProduct
0097     //
0098     // other maybe added as needs be
0099     //
0100
0101     ivhepmc = hadronizer_.getVHepMC();
0102
0103     std::vector<std::string> const& sharedResources = hadronizer_.sharedResources();
0104     for (auto const& resource : sharedResources) {
0105       usesResource(resource);
0106     }
0107
0108     if (ps.exists("ExternalDecays")) {
0109       //decayer_ = new gen::ExternalDecayDriver(ps.getParameter<ParameterSet>("ExternalDecays"));
0110       ParameterSet ps1 = ps.getParameter<ParameterSet>("ExternalDecays");
0111       decayer_ = new Decayer(ps1, consumesCollector());
0112
0113       std::vector<std::string> const& sharedResourcesDec = decayer_->sharedResources();
0114       for (auto const& resource : sharedResourcesDec) {
0115         usesResource(resource);
0116       }
0117     }
0118
0119     // This handles the case where there are no shared resources, because you
0120     // have to declare something when the SharedResources template parameter was used.
0121     if (sharedResources.empty() && (!decayer_ || decayer_->sharedResources().empty())) {
0122       usesResource(edm::uniqueSharedResourceName());
0123     }
0124
0125     if (ivhepmc == 2) {
0126       produces<edm::HepMCProduct>("unsmeared");
0127       produces<GenEventInfoProduct>();
0128     } else if (ivhepmc == 3) {
0129       produces<edm::HepMC3Product>("unsmeared");
0130       produces<GenEventInfoProduct3>();
0131     }
0132     produces<GenLumiInfoHeader, edm::Transition::BeginLuminosityBlock>();
0133     produces<GenLumiInfoProduct, edm::Transition::EndLuminosityBlock>();
0134     produces<GenRunInfoProduct, edm::Transition::EndRun>();
0135   }
0136
0137   template <class HAD, class DEC>
0138   GeneratorFilter<HAD, DEC>::~GeneratorFilter() {
0139     if (decayer_)
0140       delete decayer_;
0141   }
0142
0143   template <class HAD, class DEC>
0144   void GeneratorFilter<HAD, DEC>::preallocThreads(unsigned int iThreads) {
0145     nThreads_ = iThreads;
0146   }
0147
0148   template <class HAD, class DEC>
0149   bool GeneratorFilter<HAD, DEC>::filter(Event& ev, EventSetup const& /* es */) {
0150     RandomEngineSentry<HAD> randomEngineSentry(&hadronizer_, ev.streamID());
0151     RandomEngineSentry<DEC> randomEngineSentryDecay(decayer_, ev.streamID());
0152
0153     //added for selecting/filtering gen events, in the case of hadronizer+externalDecayer
0154
0155     bool passEvtGenSelector = false;
0156     std::unique_ptr<HepMC::GenEvent> event(nullptr);
0157     std::unique_ptr<HepMC3::GenEvent> event3(nullptr);
0158
0159     while (!passEvtGenSelector) {
0160       event.reset();
0161       event3.reset();
0162       hadronizer_.setEDMEvent(ev);
0163
0164       if (!hadronizer_.generatePartonsAndHadronize())
0165         return false;
0166
0167       //  this is "fake" stuff
0168       // in principle, decays are done as part of full event generation,
0169       // except for particles that are marked as to be kept stable
0170       // but we currently keep in it the design, because we might want
0171       // to use such feature for other applications
0172       //
0173       if (!hadronizer_.decay())
0174         return false;
0175
0176       event = hadronizer_.getGenEvent();
0177       event3 = hadronizer_.getGenEvent3();
0178       if (ivhepmc == 2 && !event.get())
0179         return false;
0180       if (ivhepmc == 3 && !event3.get())
0181         return false;
0182
0183       // The external decay driver is being added to the system,
0184       // it should be called here
0185       //
0186       if (decayer_) {
0187         if (ivhepmc == 2) {  // handle HepMC2
0188           auto t = decayer_->decay(event.get());
0189           if (t != event.get()) {
0190             event.reset(t);
0191           }
0192         }
0193         if (ivhepmc == 3) {  // handle HepMC3
0194           auto t = decayer_->decay(event3.get());
0195           if (t != event3.get()) {
0196             event3.reset(t);
0197           }
0198         }
0199       }
0200       if (ivhepmc == 2 && !event.get())
0201         return false;
0202       if (ivhepmc == 3 && !event3.get())
0203         return false;
0204
0205       passEvtGenSelector = hadronizer_.select(event.get());
0206     }
0207     // check and perform if there're any unstable particles after
0208     // running external decay packages
0209     //
0210     // fisrt of all, put back modified event tree (after external decay)
0211     //
0212     if (ivhepmc == 2)
0213       hadronizer_.resetEvent(std::move(event));
0214     else if (ivhepmc == 3)
0215       hadronizer_.resetEvent3(std::move(event3));
0216
0217     //
0218     // now run residual decays
0219     //
0220     if (!hadronizer_.residualDecay())
0221       return false;
0222
0223     hadronizer_.finalizeEvent();
0224
0225     event = hadronizer_.getGenEvent();
0226     event3 = hadronizer_.getGenEvent3();
0227     if (ivhepmc == 2) {  // HepMC
0228       if (!event.get())
0229         return false;
0230       event->set_event_number(ev.id().event());
0231
0232     } else if (ivhepmc == 3) {  // HepMC3
0233       if (!event3.get())
0234         return false;
0235       event3->set_event_number(ev.id().event());
0236     }
0237
0238     //
0239     // tutto bene - finally, form up EDM products !
0240     //
0241     if (ivhepmc == 2) {  // HepMC
0242       auto genEventInfo = hadronizer_.getGenEventInfo();
0243       if (!genEventInfo.get()) {
0244         // create GenEventInfoProduct from HepMC event in case hadronizer didn't provide one
0245         genEventInfo.reset(new GenEventInfoProduct(event.get()));
0246       }
0247
0248       ev.put(std::move(genEventInfo));
0249
0250       std::unique_ptr<HepMCProduct> bare_product(new HepMCProduct());
0251       bare_product->addHepMCData(event.release());
0252       ev.put(std::move(bare_product), "unsmeared");
0253
0254     } else if (ivhepmc == 3) {  // HepMC3
0255       auto genEventInfo3 = hadronizer_.getGenEventInfo3();
0256       if (!genEventInfo3.get()) {
0257         // create GenEventInfoProduct3 from HepMC3 event in case hadronizer didn't provide one
0258         genEventInfo3.reset(new GenEventInfoProduct3(event3.get()));
0259       }
0260
0261       ev.put(std::move(genEventInfo3));
0262
0263       std::unique_ptr<HepMC3Product> bare_product(new HepMC3Product());
0264       bare_product->addHepMCData(event3.release());
0265       ev.put(std::move(bare_product), "unsmeared");
0266     }
0267
0268     nEventsInLumiBlock_++;
0269     return true;
0270   }
0271
0272   template <class HAD, class DEC>
0273   void GeneratorFilter<HAD, DEC>::endRunProduce(Run& r, EventSetup const&) {
0274     // If relevant, record the integrated luminosity for this run
0275     // here.  To do so, we would need a standard function to invoke on
0276     // the contained hadronizer that would report the integrated
0277     // luminosity.
0278
0279     if (initialized_) {
0280       hadronizer_.statistics();
0281
0282       if (decayer_)
0283         decayer_->statistics();
0284     }
0285
0286     std::unique_ptr<GenRunInfoProduct> griproduct(new GenRunInfoProduct(hadronizer_.getGenRunInfo()));
0287     r.put(std::move(griproduct));
0288   }
0289
0290   template <class HAD, class DEC>
0291   void GeneratorFilter<HAD, DEC>::beginLuminosityBlock(LuminosityBlock const& lumi, EventSetup const& es) {}
0292
0293   template <class HAD, class DEC>
0294   void GeneratorFilter<HAD, DEC>::beginLuminosityBlockProduce(LuminosityBlock& lumi, EventSetup const& es) {
0295     nEventsInLumiBlock_ = 0;
0296     RandomEngineSentry<HAD> randomEngineSentry(&hadronizer_, lumi.index());
0297     RandomEngineSentry<DEC> randomEngineSentryDecay(decayer_, lumi.index());
0298
0299     hadronizer_.randomizeIndex(lumi, randomEngineSentry.randomEngine());
0300     hadronizer_.generateLHE(lumi, randomEngineSentry.randomEngine(), nThreads_);
0301
0302     if (!hadronizer_.readSettings(0))
0303       throw edm::Exception(errors::Configuration)
0304           << "Failed to read settings for the hadronizer " << hadronizer_.classname() << " \n";
0305
0306     if (decayer_) {
0307       decayer_->init(es);
0308       if (!hadronizer_.declareStableParticles(decayer_->operatesOnParticles()))
0309         throw edm::Exception(errors::Configuration)
0310             << "Failed to declare stable particles in hadronizer " << hadronizer_.classname() << "\n";
0311       if (!hadronizer_.declareSpecialSettings(decayer_->specialSettings()))
0312         throw edm::Exception(errors::Configuration)
0313             << "Failed to declare special settings in hadronizer " << hadronizer_.classname() << "\n";
0314     }
0315
0316     if (!hadronizer_.initializeForInternalPartons())
0317       throw edm::Exception(errors::Configuration)
0318           << "Failed to initialize hadronizer " << hadronizer_.classname() << " for internal parton generation\n";
0319
0320     std::unique_ptr<GenLumiInfoHeader> genLumiInfoHeader(hadronizer_.getGenLumiInfoHeader());
0321     lumi.put(std::move(genLumiInfoHeader));
0322     initialized_ = true;
0323   }
0324
0325   template <class HAD, class DEC>
0326   void GeneratorFilter<HAD, DEC>::endLuminosityBlock(LuminosityBlock const&, EventSetup const&) {
0327     hadronizer_.cleanLHE();
0328   }
0329
0330   template <class HAD, class DEC>
0331   void GeneratorFilter<HAD, DEC>::endLuminosityBlockProduce(LuminosityBlock& lumi, EventSetup const&) {
0332     hadronizer_.statistics();
0333     if (decayer_)
0334       decayer_->statistics();
0335
0336     GenRunInfoProduct genRunInfo = GenRunInfoProduct(hadronizer_.getGenRunInfo());
0337     std::vector<GenLumiInfoProduct::ProcessInfo> GenLumiProcess;
0338     const GenRunInfoProduct::XSec& xsec = genRunInfo.internalXSec();
0339     GenLumiInfoProduct::ProcessInfo temp;
0340     temp.setProcess(0);
0341     temp.setLheXSec(xsec.value(), xsec.error());  // Pythia gives error of -1
0342     temp.setNPassPos(nEventsInLumiBlock_);
0343     temp.setNPassNeg(0);
0344     temp.setNTotalPos(nEventsInLumiBlock_);
0345     temp.setNTotalNeg(0);
0346     temp.setTried(nEventsInLumiBlock_, nEventsInLumiBlock_, nEventsInLumiBlock_);
0347     temp.setSelected(nEventsInLumiBlock_, nEventsInLumiBlock_, nEventsInLumiBlock_);
0348     temp.setKilled(nEventsInLumiBlock_, nEventsInLumiBlock_, nEventsInLumiBlock_);
0349     temp.setAccepted(0, -1, -1);
0350     temp.setAcceptedBr(0, -1, -1);
0351     GenLumiProcess.push_back(temp);
0352
0353     std::unique_ptr<GenLumiInfoProduct> genLumiInfo(new GenLumiInfoProduct());
0354     genLumiInfo->setHEPIDWTUP(-1);
0355     genLumiInfo->setProcessInfo(GenLumiProcess);
0356
0357     lumi.put(std::move(genLumiInfo));
0358
0359     nEventsInLumiBlock_ = 0;
0360   }
0361 }  // namespace edm
0362
0363 #endif  // gen_GeneratorFilter_h