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_) {  // handle only HepMC2 for the moment
0187         auto t = decayer_->decay(event.get());
0188         if (t != event.get()) {
0189           event.reset(t);
0190         }
0191       }
0192       if (ivhepmc == 2 && !event.get())
0193         return false;
0194       if (ivhepmc == 3 && !event3.get())
0195         return false;
0196
0197       passEvtGenSelector = hadronizer_.select(event.get());
0198     }
0199     // check and perform if there're any unstable particles after
0200     // running external decay packages
0201     //
0202     // fisrt of all, put back modified event tree (after external decay)
0203     //
0204     if (ivhepmc == 2)
0205       hadronizer_.resetEvent(std::move(event));
0206     else if (ivhepmc == 3)
0207       hadronizer_.resetEvent3(std::move(event3));
0208
0209     //
0210     // now run residual decays
0211     //
0212     if (!hadronizer_.residualDecay())
0213       return false;
0214
0215     hadronizer_.finalizeEvent();
0216
0217     event = hadronizer_.getGenEvent();
0218     event3 = hadronizer_.getGenEvent3();
0219     if (ivhepmc == 2) {  // HepMC
0220       if (!event.get())
0221         return false;
0222       event->set_event_number(ev.id().event());
0223
0224     } else if (ivhepmc == 3) {  // HepMC3
0225       if (!event3.get())
0226         return false;
0227       event3->set_event_number(ev.id().event());
0228     }
0229
0230     //
0231     // tutto bene - finally, form up EDM products !
0232     //
0233     if (ivhepmc == 2) {  // HepMC
0234       auto genEventInfo = hadronizer_.getGenEventInfo();
0235       if (!genEventInfo.get()) {
0236         // create GenEventInfoProduct from HepMC event in case hadronizer didn't provide one
0237         genEventInfo.reset(new GenEventInfoProduct(event.get()));
0238       }
0239
0240       ev.put(std::move(genEventInfo));
0241
0242       std::unique_ptr<HepMCProduct> bare_product(new HepMCProduct());
0243       bare_product->addHepMCData(event.release());
0244       ev.put(std::move(bare_product), "unsmeared");
0245
0246     } else if (ivhepmc == 3) {  // HepMC3
0247       auto genEventInfo3 = hadronizer_.getGenEventInfo3();
0248       if (!genEventInfo3.get()) {
0249         // create GenEventInfoProduct3 from HepMC3 event in case hadronizer didn't provide one
0250         genEventInfo3.reset(new GenEventInfoProduct3(event3.get()));
0251       }
0252
0253       ev.put(std::move(genEventInfo3));
0254
0255       std::unique_ptr<HepMC3Product> bare_product(new HepMC3Product());
0256       bare_product->addHepMCData(event3.release());
0257       ev.put(std::move(bare_product), "unsmeared");
0258     }
0259
0260     nEventsInLumiBlock_++;
0261     return true;
0262   }
0263
0264   template <class HAD, class DEC>
0265   void GeneratorFilter<HAD, DEC>::endRunProduce(Run& r, EventSetup const&) {
0266     // If relevant, record the integrated luminosity for this run
0267     // here.  To do so, we would need a standard function to invoke on
0268     // the contained hadronizer that would report the integrated
0269     // luminosity.
0270
0271     if (initialized_) {
0272       hadronizer_.statistics();
0273
0274       if (decayer_)
0275         decayer_->statistics();
0276     }
0277
0278     std::unique_ptr<GenRunInfoProduct> griproduct(new GenRunInfoProduct(hadronizer_.getGenRunInfo()));
0279     r.put(std::move(griproduct));
0280   }
0281
0282   template <class HAD, class DEC>
0283   void GeneratorFilter<HAD, DEC>::beginLuminosityBlock(LuminosityBlock const& lumi, EventSetup const& es) {}
0284
0285   template <class HAD, class DEC>
0286   void GeneratorFilter<HAD, DEC>::beginLuminosityBlockProduce(LuminosityBlock& lumi, EventSetup const& es) {
0287     nEventsInLumiBlock_ = 0;
0288     RandomEngineSentry<HAD> randomEngineSentry(&hadronizer_, lumi.index());
0289     RandomEngineSentry<DEC> randomEngineSentryDecay(decayer_, lumi.index());
0290
0291     hadronizer_.randomizeIndex(lumi, randomEngineSentry.randomEngine());
0292     hadronizer_.generateLHE(lumi, randomEngineSentry.randomEngine(), nThreads_);
0293
0294     if (!hadronizer_.readSettings(0))
0295       throw edm::Exception(errors::Configuration)
0296           << "Failed to read settings for the hadronizer " << hadronizer_.classname() << " \n";
0297
0298     if (decayer_) {
0299       decayer_->init(es);
0300       if (!hadronizer_.declareStableParticles(decayer_->operatesOnParticles()))
0301         throw edm::Exception(errors::Configuration)
0302             << "Failed to declare stable particles in hadronizer " << hadronizer_.classname() << "\n";
0303       if (!hadronizer_.declareSpecialSettings(decayer_->specialSettings()))
0304         throw edm::Exception(errors::Configuration)
0305             << "Failed to declare special settings in hadronizer " << hadronizer_.classname() << "\n";
0306     }
0307
0308     if (!hadronizer_.initializeForInternalPartons())
0309       throw edm::Exception(errors::Configuration)
0310           << "Failed to initialize hadronizer " << hadronizer_.classname() << " for internal parton generation\n";
0311
0312     std::unique_ptr<GenLumiInfoHeader> genLumiInfoHeader(hadronizer_.getGenLumiInfoHeader());
0313     lumi.put(std::move(genLumiInfoHeader));
0314     initialized_ = true;
0315   }
0316
0317   template <class HAD, class DEC>
0318   void GeneratorFilter<HAD, DEC>::endLuminosityBlock(LuminosityBlock const&, EventSetup const&) {
0319     hadronizer_.cleanLHE();
0320   }
0321
0322   template <class HAD, class DEC>
0323   void GeneratorFilter<HAD, DEC>::endLuminosityBlockProduce(LuminosityBlock& lumi, EventSetup const&) {
0324     hadronizer_.statistics();
0325     if (decayer_)
0326       decayer_->statistics();
0327
0328     GenRunInfoProduct genRunInfo = GenRunInfoProduct(hadronizer_.getGenRunInfo());
0329     std::vector<GenLumiInfoProduct::ProcessInfo> GenLumiProcess;
0330     const GenRunInfoProduct::XSec& xsec = genRunInfo.internalXSec();
0331     GenLumiInfoProduct::ProcessInfo temp;
0332     temp.setProcess(0);
0333     temp.setLheXSec(xsec.value(), xsec.error());  // Pythia gives error of -1
0334     temp.setNPassPos(nEventsInLumiBlock_);
0335     temp.setNPassNeg(0);
0336     temp.setNTotalPos(nEventsInLumiBlock_);
0337     temp.setNTotalNeg(0);
0338     temp.setTried(nEventsInLumiBlock_, nEventsInLumiBlock_, nEventsInLumiBlock_);
0339     temp.setSelected(nEventsInLumiBlock_, nEventsInLumiBlock_, nEventsInLumiBlock_);
0340     temp.setKilled(nEventsInLumiBlock_, nEventsInLumiBlock_, nEventsInLumiBlock_);
0341     temp.setAccepted(0, -1, -1);
0342     temp.setAcceptedBr(0, -1, -1);
0343     GenLumiProcess.push_back(temp);
0344
0345     std::unique_ptr<GenLumiInfoProduct> genLumiInfo(new GenLumiInfoProduct());
0346     genLumiInfo->setHEPIDWTUP(-1);
0347     genLumiInfo->setProcessInfo(GenLumiProcess);
0348
0349     lumi.put(std::move(genLumiInfo));
0350
0351     nEventsInLumiBlock_ = 0;
0352   }
0353 }  // namespace edm
0354
0355 #endif  // gen_GeneratorFilter_h