Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-06-11 03:13:11

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