Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-02-28 23:29:41

0001 #include <iostream>
0002 #include <memory>
0003 
0004 #include <cstdint>
0005 #include <memory>
0006 #include <sstream>
0007 #include <string>
0008 #include <vector>
0009 
0010 #include "HepMC/GenEvent.h"
0011 #include "HepMC/GenParticle.h"
0012 
0013 #include "Pythia8/Pythia.h"
0014 #include "Pythia8Plugins/HepMC2.h"
0015 
0016 using namespace Pythia8;
0017 
0018 #include "GeneratorInterface/Pythia8Interface/interface/Py8InterfaceBase.h"
0019 
0020 #include "ReweightUserHooks.h"
0021 #include "GeneratorInterface/Pythia8Interface/interface/CustomHook.h"
0022 #include "TopRecoilHook.h"
0023 
0024 // PS matchning prototype
0025 //
0026 #include "GeneratorInterface/Pythia8Interface/plugins/JetMatchingHook.h"
0027 #include "Pythia8Plugins/JetMatching.h"
0028 #include "Pythia8Plugins/aMCatNLOHooks.h"
0029 
0030 // Emission Veto Hooks
0031 //
0032 #include "Pythia8Plugins/PowhegHooks.h"
0033 #include "GeneratorInterface/Pythia8Interface/plugins/EmissionVetoHook1.h"
0034 
0035 // Resonance scale hook
0036 #include "GeneratorInterface/Pythia8Interface/plugins/PowhegResHook.h"
0037 #include "GeneratorInterface/Pythia8Interface/plugins/PowhegHooksBB4L.h"
0038 
0039 //biased tau decayer
0040 #include "GeneratorInterface/Pythia8Interface/interface/BiasedTauDecayer.h"
0041 
0042 //decay filter hook
0043 #include "GeneratorInterface/Pythia8Interface/interface/ResonanceDecayFilterHook.h"
0044 
0045 //decay filter hook
0046 #include "GeneratorInterface/Pythia8Interface/interface/PTFilterHook.h"
0047 
0048 // EvtGen plugin
0049 //
0050 #include "Pythia8Plugins/EvtGen.h"
0051 
0052 #include "FWCore/Concurrency/interface/SharedResourceNames.h"
0053 #include "FWCore/ServiceRegistry/interface/Service.h"
0054 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0055 #include "FWCore/Utilities/interface/RandomNumberGenerator.h"
0056 #include "FWCore/ParameterSet/interface/FileInPath.h"
0057 
0058 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0059 #include "SimDataFormats/GeneratorProducts/interface/GenRunInfoProduct.h"
0060 
0061 #include "GeneratorInterface/Core/interface/GeneratorFilter.h"
0062 #include "GeneratorInterface/Core/interface/HadronizerFilter.h"
0063 #include "GeneratorInterface/Core/interface/ConcurrentGeneratorFilter.h"
0064 #include "GeneratorInterface/Core/interface/ConcurrentHadronizerFilter.h"
0065 
0066 #include "GeneratorInterface/Pythia8Interface/plugins/LHAupLesHouches.h"
0067 
0068 #include "HepPID/ParticleIDTranslations.hh"
0069 
0070 #include "GeneratorInterface/ExternalDecays/interface/ExternalDecayDriver.h"
0071 #include "GeneratorInterface/ExternalDecays/interface/ConcurrentExternalDecayDriver.h"
0072 
0073 namespace CLHEP {
0074   class HepRandomEngine;
0075 }
0076 
0077 using namespace gen;
0078 
0079 //Insert class for use w/ PDFPtr for proton-photon flux
0080 //parameters hardcoded according to main70.cc in PYTHIA8 v3.10
0081 class Nucleus2gamma2 : public Pythia8::PDF {
0082 public:
0083   // Constructor.
0084   Nucleus2gamma2(int idBeamIn) : Pythia8::PDF(idBeamIn) {}
0085 
0086   // Update the photon flux.
0087   void xfUpdate(int idBeamIn, double x, double) override {
0088     // lead
0089     double radius = 0;  // radius in [fm]
0090     double z = 0;
0091     if (idBeamIn == 1000822080) {
0092       radius = 6.636;
0093       z = 82;
0094     }
0095     // oxygen
0096     else if (idBeamIn == 80160) {
0097       radius = 3.02;
0098       z = 8;
0099     }
0100 
0101     // Minimum impact parameter (~2*radius) [fm].
0102     double bmin = 2 * radius;
0103 
0104     // Per-nucleon mass for lead.
0105     double m2 = pow2(0.9314);
0106     double alphaEM = 0.007297353080;
0107     double hbarc = 0.197;
0108     double xi = x * sqrt(m2) * bmin / hbarc;
0109     double bK0 = besselK0(xi);
0110     double bK1 = besselK1(xi);
0111     double intB = xi * bK1 * bK0 - 0.5 * pow2(xi) * (pow2(bK1) - pow2(bK0));
0112     xgamma = 2. * alphaEM * pow2(z) / M_PI * intB;
0113   }
0114 };
0115 
0116 class Pythia8Hadronizer : public Py8InterfaceBase {
0117 public:
0118   Pythia8Hadronizer(const edm::ParameterSet &params);
0119   ~Pythia8Hadronizer() override;
0120 
0121   bool initializeForInternalPartons() override;
0122   bool initializeForExternalPartons();
0123 
0124   bool generatePartonsAndHadronize() override;
0125   bool hadronize();
0126 
0127   virtual bool residualDecay();
0128 
0129   void finalizeEvent() override;
0130 
0131   void statistics() override;
0132 
0133   const char *classname() const override { return "Pythia8Hadronizer"; }
0134 
0135   std::unique_ptr<GenLumiInfoHeader> getGenLumiInfoHeader() const override;
0136 
0137 private:
0138   void doSetRandomEngine(CLHEP::HepRandomEngine *v) override { p8SetRandomEngine(v); }
0139   std::vector<std::string> const &doSharedResources() const override { return p8SharedResources; }
0140 
0141   /// Center-of-Mass energy
0142   double comEnergy;
0143 
0144   std::string LHEInputFileName;
0145   std::shared_ptr<LHAupLesHouches> lhaUP;
0146 
0147   enum { PP, PPbar, ElectronPositron };
0148   int fInitialState;  // pp, ppbar, or e-e+
0149 
0150   double fBeam1PZ;
0151   double fBeam2PZ;
0152 
0153   //PDFPtr for the photonFlux
0154   //Following main70.cc example in PYTHIA8 v3.10
0155   bool doProtonPhotonFlux = false;
0156   Pythia8::PDFPtr photonFlux = nullptr;
0157 
0158   //helper class to allow multiple user hooks simultaneously
0159   std::shared_ptr<UserHooksVector> fUserHooksVector;
0160   bool UserHooksSet;
0161 
0162   // Reweight user hooks
0163   //
0164   std::shared_ptr<UserHooks> fReweightUserHook;
0165   std::shared_ptr<UserHooks> fReweightEmpUserHook;
0166   std::shared_ptr<UserHooks> fReweightRapUserHook;
0167   std::shared_ptr<UserHooks> fReweightPtHatRapUserHook;
0168 
0169   // PS matching prototype
0170   //
0171   std::shared_ptr<JetMatchingHook> fJetMatchingHook;
0172   std::shared_ptr<Pythia8::JetMatchingMadgraph> fJetMatchingPy8InternalHook;
0173   std::shared_ptr<Pythia8::amcnlo_unitarised_interface> fMergingHook;
0174 
0175   // Emission Veto Hooks
0176   //
0177   std::shared_ptr<PowhegHooks> fEmissionVetoHook;
0178   std::shared_ptr<EmissionVetoHook1> fEmissionVetoHook1;
0179 
0180   // Resonance scale hook
0181   std::shared_ptr<PowhegResHook> fPowhegResHook;
0182   std::shared_ptr<PowhegHooksBB4L> fPowhegHooksBB4L;
0183 
0184   // biased tau decayer
0185   std::shared_ptr<BiasedTauDecayer> fBiasedTauDecayer;
0186 
0187   //resonance decay filter hook
0188   std::shared_ptr<ResonanceDecayFilterHook> fResonanceDecayFilterHook;
0189 
0190   //PT filter hook
0191   std::shared_ptr<PTFilterHook> fPTFilterHook;
0192 
0193   //Generic customized hooks vector
0194   std::shared_ptr<UserHooksVector> fCustomHooksVector;
0195 
0196   //RecoilToTop userhook
0197   std::shared_ptr<TopRecoilHook> fTopRecoilHook;
0198 
0199   int EV1_nFinal;
0200   bool EV1_vetoOn;
0201   int EV1_maxVetoCount;
0202   int EV1_pThardMode;
0203   int EV1_pTempMode;
0204   int EV1_emittedMode;
0205   int EV1_pTdefMode;
0206   bool EV1_MPIvetoOn;
0207   int EV1_QEDvetoMode;
0208   int EV1_nFinalMode;
0209 
0210   static const std::vector<std::string> p8SharedResources;
0211 
0212   vector<float> DJR;
0213   int nME;
0214   int nMEFiltered;
0215 
0216   int nISRveto;
0217   int nFSRveto;
0218 };
0219 
0220 const std::vector<std::string> Pythia8Hadronizer::p8SharedResources = {edm::SharedResourceNames::kPythia8};
0221 
0222 Pythia8Hadronizer::Pythia8Hadronizer(const edm::ParameterSet &params)
0223     : Py8InterfaceBase(params),
0224       comEnergy(params.getParameter<double>("comEnergy")),
0225       LHEInputFileName(params.getUntrackedParameter<std::string>("LHEInputFileName", "")),
0226       fInitialState(PP),
0227       doProtonPhotonFlux(params.getUntrackedParameter<bool>("doProtonPhotonFlux", false)),
0228       UserHooksSet(false),
0229       nME(-1),
0230       nMEFiltered(-1),
0231       nISRveto(0),
0232       nFSRveto(0) {
0233   ivhepmc = 2;
0234   // J.Y.: the following 3 parameters are hacked "for a reason"
0235   //
0236   if (params.exists("PPbarInitialState")) {
0237     if (fInitialState == PP) {
0238       fInitialState = PPbar;
0239       edm::LogImportant("GeneratorInterface|Pythia8Interface")
0240           << "Pythia8 will be initialized for PROTON-ANTIPROTON INITIAL STATE. "
0241           << "This is a user-request change from the DEFAULT PROTON-PROTON initial state.";
0242     } else {
0243       // probably need to throw on attempt to override ?
0244     }
0245   } else if (params.exists("ElectronPositronInitialState")) {
0246     if (fInitialState == PP) {
0247       fInitialState = ElectronPositron;
0248       edm::LogInfo("GeneratorInterface|Pythia8Interface")
0249           << "Pythia8 will be initialized for ELECTRON-POSITRON INITIAL STATE. "
0250           << "This is a user-request change from the DEFAULT PROTON-PROTON initial state.";
0251     } else {
0252       // probably need to throw on attempt to override ?
0253     }
0254   } else if (params.exists("ElectronProtonInitialState") || params.exists("PositronProtonInitialState")) {
0255     // throw on unknown initial state !
0256     throw edm::Exception(edm::errors::Configuration, "Pythia8Interface")
0257         << " UNKNOWN INITIAL STATE. \n The allowed initial states are: PP, PPbar, ElectronPositron \n";
0258   }
0259 
0260   // avoid filling weights twice (from v8.30x)
0261   toHepMC.set_store_weights(false);
0262 
0263   // Reweight user hook
0264   //
0265   if (params.exists("reweightGen")) {
0266     edm::LogInfo("Pythia8Interface") << "Start setup for reweightGen";
0267     edm::ParameterSet rgParams = params.getParameter<edm::ParameterSet>("reweightGen");
0268     fReweightUserHook.reset(
0269         new PtHatReweightUserHook(rgParams.getParameter<double>("pTRef"), rgParams.getParameter<double>("power")));
0270     edm::LogInfo("Pythia8Interface") << "End setup for reweightGen";
0271   }
0272   if (params.exists("reweightGenEmp")) {
0273     edm::LogInfo("Pythia8Interface") << "Start setup for reweightGenEmp";
0274     edm::ParameterSet rgeParams = params.getParameter<edm::ParameterSet>("reweightGenEmp");
0275 
0276     std::string tuneName = "";
0277     if (rgeParams.exists("tune"))
0278       tuneName = rgeParams.getParameter<std::string>("tune");
0279     fReweightEmpUserHook.reset(new PtHatEmpReweightUserHook(tuneName));
0280     edm::LogInfo("Pythia8Interface") << "End setup for reweightGenEmp";
0281   }
0282   if (params.exists("reweightGenRap")) {
0283     edm::LogInfo("Pythia8Interface") << "Start setup for reweightGenRap";
0284     edm::ParameterSet rgrParams = params.getParameter<edm::ParameterSet>("reweightGenRap");
0285     fReweightRapUserHook.reset(new RapReweightUserHook(rgrParams.getParameter<std::string>("yLabSigmaFunc"),
0286                                                        rgrParams.getParameter<double>("yLabPower"),
0287                                                        rgrParams.getParameter<std::string>("yCMSigmaFunc"),
0288                                                        rgrParams.getParameter<double>("yCMPower"),
0289                                                        rgrParams.getParameter<double>("pTHatMin"),
0290                                                        rgrParams.getParameter<double>("pTHatMax")));
0291     edm::LogInfo("Pythia8Interface") << "End setup for reweightGenRap";
0292   }
0293   if (params.exists("reweightGenPtHatRap")) {
0294     edm::LogInfo("Pythia8Interface") << "Start setup for reweightGenPtHatRap";
0295     edm::ParameterSet rgrParams = params.getParameter<edm::ParameterSet>("reweightGenPtHatRap");
0296     fReweightPtHatRapUserHook.reset(new PtHatRapReweightUserHook(rgrParams.getParameter<std::string>("yLabSigmaFunc"),
0297                                                                  rgrParams.getParameter<double>("yLabPower"),
0298                                                                  rgrParams.getParameter<std::string>("yCMSigmaFunc"),
0299                                                                  rgrParams.getParameter<double>("yCMPower"),
0300                                                                  rgrParams.getParameter<double>("pTHatMin"),
0301                                                                  rgrParams.getParameter<double>("pTHatMax")));
0302     edm::LogInfo("Pythia8Interface") << "End setup for reweightGenPtHatRap";
0303   }
0304 
0305   if (params.exists("useUserHook"))
0306     throw edm::Exception(edm::errors::Configuration, "Pythia8Interface")
0307         << " Obsolete parameter: useUserHook \n Please use the actual one instead \n";
0308 
0309   // PS matching prototype
0310   //
0311   if (params.exists("jetMatching")) {
0312     edm::ParameterSet jmParams = params.getUntrackedParameter<edm::ParameterSet>("jetMatching");
0313     std::string scheme = jmParams.getParameter<std::string>("scheme");
0314     if (scheme == "Madgraph" || scheme == "MadgraphFastJet") {
0315       fJetMatchingHook.reset(new JetMatchingHook(jmParams, &fMasterGen->info));
0316     }
0317   }
0318 
0319   // Pythia8Interface emission veto
0320   //
0321   if (params.exists("emissionVeto1")) {
0322     EV1_nFinal = -1;
0323     if (params.exists("EV1_nFinal"))
0324       EV1_nFinal = params.getParameter<int>("EV1_nFinal");
0325     EV1_vetoOn = true;
0326     if (params.exists("EV1_vetoOn"))
0327       EV1_vetoOn = params.getParameter<bool>("EV1_vetoOn");
0328     EV1_maxVetoCount = 10;
0329     if (params.exists("EV1_maxVetoCount"))
0330       EV1_maxVetoCount = params.getParameter<int>("EV1_maxVetoCount");
0331     EV1_pThardMode = 1;
0332     if (params.exists("EV1_pThardMode"))
0333       EV1_pThardMode = params.getParameter<int>("EV1_pThardMode");
0334     EV1_pTempMode = 0;
0335     if (params.exists("EV1_pTempMode"))
0336       EV1_pTempMode = params.getParameter<int>("EV1_pTempMode");
0337     if (EV1_pTempMode > 2 || EV1_pTempMode < 0)
0338       throw edm::Exception(edm::errors::Configuration, "Pythia8Interface") << " Wrong value for EV1_pTempMode code\n";
0339     EV1_emittedMode = 0;
0340     if (params.exists("EV1_emittedMode"))
0341       EV1_emittedMode = params.getParameter<int>("EV1_emittedMode");
0342     EV1_pTdefMode = 1;
0343     if (params.exists("EV1_pTdefMode"))
0344       EV1_pTdefMode = params.getParameter<int>("EV1_pTdefMode");
0345     EV1_MPIvetoOn = false;
0346     if (params.exists("EV1_MPIvetoOn"))
0347       EV1_MPIvetoOn = params.getParameter<bool>("EV1_MPIvetoOn");
0348     EV1_QEDvetoMode = 0;
0349     if (params.exists("EV1_QEDvetoMode"))
0350       EV1_QEDvetoMode = params.getParameter<int>("EV1_QEDvetoMode");
0351     EV1_nFinalMode = 0;
0352     if (params.exists("EV1_nFinalMode"))
0353       EV1_nFinalMode = params.getParameter<int>("EV1_nFinalMode");
0354     fEmissionVetoHook1.reset(new EmissionVetoHook1(EV1_nFinal,
0355                                                    EV1_vetoOn,
0356                                                    EV1_maxVetoCount,
0357                                                    EV1_pThardMode,
0358                                                    EV1_pTempMode,
0359                                                    EV1_emittedMode,
0360                                                    EV1_pTdefMode,
0361                                                    EV1_MPIvetoOn,
0362                                                    EV1_QEDvetoMode,
0363                                                    EV1_nFinalMode,
0364                                                    0));
0365   }
0366 
0367   if (params.exists("UserCustomization")) {
0368     fCustomHooksVector = std::make_shared<UserHooksVector>();
0369     const std::vector<edm::ParameterSet> userParams =
0370         params.getParameter<std::vector<edm::ParameterSet>>("UserCustomization");
0371     for (const auto &pluginParams : userParams) {
0372       (fCustomHooksVector->hooks)
0373           .push_back(
0374               CustomHookFactory::get()->create(pluginParams.getParameter<std::string>("pluginName"), pluginParams));
0375     }
0376   }
0377 
0378   if (params.exists("VinciaPlugin")) {
0379     throw edm::Exception(edm::errors::Configuration, "Pythia8Interface")
0380         << " Obsolete parameter: VinciaPlugin \n Please use the parameter PartonShowers:model instead \n";
0381   }
0382   if (params.exists("DirePlugin")) {
0383     throw edm::Exception(edm::errors::Configuration, "Pythia8Interface")
0384         << " Obsolete parameter: DirePlugin \n Please use the parameter PartonShowers:model instead \n";
0385   }
0386 }
0387 
0388 Pythia8Hadronizer::~Pythia8Hadronizer() {}
0389 
0390 bool Pythia8Hadronizer::initializeForInternalPartons() {
0391   bool status = false, status1 = false;
0392 
0393   if (lheFile_.empty()) {
0394     if (fInitialState == PP)  // default
0395     {
0396       fMasterGen->settings.mode("Beams:idA", 2212);
0397       fMasterGen->settings.mode("Beams:idB", 2212);
0398     } else if (fInitialState == PPbar) {
0399       fMasterGen->settings.mode("Beams:idA", 2212);
0400       fMasterGen->settings.mode("Beams:idB", -2212);
0401     } else if (fInitialState == ElectronPositron) {
0402       fMasterGen->settings.mode("Beams:idA", 11);
0403       fMasterGen->settings.mode("Beams:idB", -11);
0404     } else {
0405       // throw on unknown initial state !
0406       throw edm::Exception(edm::errors::Configuration, "Pythia8Interface")
0407           << " UNKNOWN INITIAL STATE. \n The allowed initial states are: PP, PPbar, ElectronPositron \n";
0408     }
0409     fMasterGen->settings.parm("Beams:eCM", comEnergy);
0410   } else {
0411     fMasterGen->settings.mode("Beams:frameType", 4);
0412     fMasterGen->settings.word("Beams:LHEF", lheFile_);
0413   }
0414 
0415   if (doProtonPhotonFlux) {
0416     photonFlux = make_shared<Nucleus2gamma2>(2212);
0417     fMasterGen->setPhotonFluxPtr(photonFlux, nullptr);
0418   }
0419 
0420   if (!fUserHooksVector.get())
0421     fUserHooksVector.reset(new UserHooksVector);
0422   (fUserHooksVector->hooks).clear();
0423 
0424   if (fReweightUserHook.get())
0425     (fUserHooksVector->hooks).push_back(fReweightUserHook);
0426   if (fReweightEmpUserHook.get())
0427     (fUserHooksVector->hooks).push_back(fReweightEmpUserHook);
0428   if (fReweightRapUserHook.get())
0429     (fUserHooksVector->hooks).push_back(fReweightRapUserHook);
0430   if (fReweightPtHatRapUserHook.get())
0431     (fUserHooksVector->hooks).push_back(fReweightPtHatRapUserHook);
0432   if (fJetMatchingHook.get())
0433     (fUserHooksVector->hooks).push_back(fJetMatchingHook);
0434   if (fEmissionVetoHook1.get()) {
0435     edm::LogInfo("Pythia8Interface") << "Turning on Emission Veto Hook 1 from CMSSW Pythia8Interface";
0436     (fUserHooksVector->hooks).push_back(fEmissionVetoHook1);
0437   }
0438 
0439   if (fMasterGen->settings.mode("POWHEG:veto") > 0 || fMasterGen->settings.mode("POWHEG:MPIveto") > 0) {
0440     if (fJetMatchingHook.get() || fEmissionVetoHook1.get())
0441       throw edm::Exception(edm::errors::Configuration, "Pythia8Interface")
0442           << " Attempt to turn on PowhegHooks by pythia8 settings but there are incompatible hooks on \n Incompatible "
0443              "are : jetMatching, emissionVeto1 \n";
0444 
0445     if (!fEmissionVetoHook.get())
0446       fEmissionVetoHook.reset(new PowhegHooks());
0447 
0448     edm::LogInfo("Pythia8Interface") << "Turning on Emission Veto Hook from pythia8 code";
0449     (fUserHooksVector->hooks).push_back(fEmissionVetoHook);
0450   }
0451 
0452   bool PowhegRes = fMasterGen->settings.flag("POWHEGres:calcScales");
0453   if (PowhegRes) {
0454     edm::LogInfo("Pythia8Interface") << "Turning on resonance scale setting from CMSSW Pythia8Interface";
0455     if (!fPowhegResHook.get())
0456       fPowhegResHook.reset(new PowhegResHook());
0457     (fUserHooksVector->hooks).push_back(fPowhegResHook);
0458   }
0459 
0460   bool PowhegBB4L = fMasterGen->settings.flag("POWHEG:bb4l");
0461   if (PowhegBB4L) {
0462     edm::LogInfo("Pythia8Interface") << "Turning on BB4l hook from CMSSW Pythia8Interface";
0463     if (!fPowhegHooksBB4L.get())
0464       fPowhegHooksBB4L.reset(new PowhegHooksBB4L());
0465     (fUserHooksVector->hooks).push_back(fPowhegHooksBB4L);
0466   }
0467 
0468   bool TopRecoilHook1 = fMasterGen->settings.flag("TopRecoilHook:doTopRecoilIn");
0469   if (TopRecoilHook1) {
0470     edm::LogInfo("Pythia8Interface") << "Turning on RecoilToTop hook from Pythia8Interface";
0471     if (!fTopRecoilHook.get())
0472       fTopRecoilHook.reset(new TopRecoilHook());
0473     (fUserHooksVector->hooks).push_back(fTopRecoilHook);
0474   }
0475 
0476   //adapted from main89.cc in pythia8 examples
0477   bool internalMatching = fMasterGen->settings.flag("JetMatching:merge");
0478   bool internalMerging = !(fMasterGen->settings.word("Merging:Process") == "void");
0479 
0480   if (internalMatching && internalMerging) {
0481     throw edm::Exception(edm::errors::Configuration, "Pythia8Interface")
0482         << " Only one jet matching/merging scheme can be used at a time. \n";
0483   }
0484 
0485   if (internalMatching) {
0486     if (!fJetMatchingPy8InternalHook.get())
0487       fJetMatchingPy8InternalHook.reset(new Pythia8::JetMatchingMadgraph);
0488     (fUserHooksVector->hooks).push_back(fJetMatchingPy8InternalHook);
0489   }
0490 
0491   if (internalMerging) {
0492     int scheme = (fMasterGen->settings.flag("Merging:doUMEPSTree") || fMasterGen->settings.flag("Merging:doUMEPSSubt"))
0493                      ? 1
0494                      : ((fMasterGen->settings.flag("Merging:doUNLOPSTree") ||
0495                          fMasterGen->settings.flag("Merging:doUNLOPSSubt") ||
0496                          fMasterGen->settings.flag("Merging:doUNLOPSLoop") ||
0497                          fMasterGen->settings.flag("Merging:doUNLOPSSubtNLO"))
0498                             ? 2
0499                             : 0);
0500     if (!fMergingHook.get())
0501       fMergingHook.reset(new Pythia8::amcnlo_unitarised_interface(scheme));
0502     (fUserHooksVector->hooks).push_back(fMergingHook);
0503   }
0504 
0505   bool biasedTauDecayer = fMasterGen->settings.flag("BiasedTauDecayer:filter");
0506   if (biasedTauDecayer) {
0507     if (!fBiasedTauDecayer.get()) {
0508       Pythia8::Info localInfo = fMasterGen->info;
0509       fBiasedTauDecayer.reset(new BiasedTauDecayer(&localInfo, &(fMasterGen->settings)));
0510     }
0511     std::vector<int> handledParticles;
0512     handledParticles.push_back(15);
0513     fMasterGen->setDecayPtr(fBiasedTauDecayer, handledParticles);
0514   }
0515 
0516   bool resonanceDecayFilter = fMasterGen->settings.flag("ResonanceDecayFilter:filter");
0517   if (resonanceDecayFilter) {
0518     fResonanceDecayFilterHook.reset(new ResonanceDecayFilterHook);
0519     (fUserHooksVector->hooks).push_back(fResonanceDecayFilterHook);
0520   }
0521 
0522   bool PTFilter = fMasterGen->settings.flag("PTFilter:filter");
0523   if (PTFilter) {
0524     fPTFilterHook.reset(new PTFilterHook);
0525     (fUserHooksVector->hooks).push_back(fPTFilterHook);
0526   }
0527 
0528   if (!(fUserHooksVector->hooks).empty() && !UserHooksSet) {
0529     for (auto &fUserHook : fUserHooksVector->hooks) {
0530       fMasterGen->addUserHooksPtr(fUserHook);
0531     }
0532     UserHooksSet = true;
0533   }
0534 
0535   if (fCustomHooksVector.get()) {
0536     edm::LogInfo("Pythia8Interface") << "Adding customized user hooks";
0537     for (const auto &fUserHook : fCustomHooksVector->hooks) {
0538       fMasterGen->addUserHooksPtr(fUserHook);
0539     }
0540   }
0541 
0542   edm::LogInfo("Pythia8Interface") << "Initializing MasterGen";
0543   status = fMasterGen->init();
0544 
0545   //clean up temp file
0546   if (!slhafile_.empty()) {
0547     std::remove(slhafile_.c_str());
0548   }
0549 
0550   if (pythiaPylistVerbosity > 10) {
0551     if (pythiaPylistVerbosity == 11 || pythiaPylistVerbosity == 13)
0552       fMasterGen->settings.listAll();
0553     if (pythiaPylistVerbosity == 12 || pythiaPylistVerbosity == 13)
0554       fMasterGen->particleData.listAll();
0555   }
0556 
0557   // init decayer
0558   fDecayer->settings.flag("ProcessLevel:all", false);  // trick
0559   fDecayer->settings.flag("ProcessLevel:resonanceDecays", true);
0560   edm::LogInfo("Pythia8Interface") << "Initializing Decayer";
0561   status1 = fDecayer->init();
0562 
0563   if (useEvtGen) {
0564     edm::LogInfo("Pythia8Hadronizer") << "Creating and initializing pythia8 EvtGen plugin";
0565     if (!evtgenDecays.get()) {
0566       evtgenDecays.reset(new EvtGenDecays(fMasterGen.get(), evtgenDecFile, evtgenPdlFile));
0567       for (unsigned int i = 0; i < evtgenUserFiles.size(); i++)
0568         evtgenDecays->readDecayFile(evtgenUserFiles.at(i));
0569     }
0570   }
0571 
0572   return (status && status1);
0573 }
0574 
0575 bool Pythia8Hadronizer::initializeForExternalPartons() {
0576   edm::LogInfo("Pythia8Interface") << "Initializing for external partons";
0577 
0578   bool status = false, status1 = false;
0579 
0580   if (!fUserHooksVector.get())
0581     fUserHooksVector.reset(new UserHooksVector);
0582   (fUserHooksVector->hooks).clear();
0583 
0584   if (fReweightUserHook.get())
0585     (fUserHooksVector->hooks).push_back(fReweightUserHook);
0586   if (fReweightEmpUserHook.get())
0587     (fUserHooksVector->hooks).push_back(fReweightEmpUserHook);
0588   if (fReweightRapUserHook.get())
0589     (fUserHooksVector->hooks).push_back(fReweightRapUserHook);
0590   if (fReweightPtHatRapUserHook.get())
0591     (fUserHooksVector->hooks).push_back(fReweightPtHatRapUserHook);
0592   if (fJetMatchingHook.get())
0593     (fUserHooksVector->hooks).push_back(fJetMatchingHook);
0594   if (fEmissionVetoHook1.get()) {
0595     edm::LogInfo("Pythia8Interface") << "Turning on Emission Veto Hook 1 from CMSSW Pythia8Interface";
0596     (fUserHooksVector->hooks).push_back(fEmissionVetoHook1);
0597   }
0598 
0599   if (fCustomHooksVector.get()) {
0600     edm::LogInfo("Pythia8Interface") << "Adding customized user hook";
0601     for (const auto &fUserHook : fCustomHooksVector->hooks) {
0602       (fUserHooksVector->hooks).push_back(fUserHook);
0603     }
0604   }
0605 
0606   if (fMasterGen->settings.mode("POWHEG:veto") > 0 || fMasterGen->settings.mode("POWHEG:MPIveto") > 0) {
0607     if (fJetMatchingHook.get() || fEmissionVetoHook1.get())
0608       throw edm::Exception(edm::errors::Configuration, "Pythia8Interface")
0609           << " Attempt to turn on PowhegHooks by pythia8 settings but there are incompatible hooks on \n Incompatible "
0610              "are : jetMatching, emissionVeto1 \n";
0611 
0612     if (!fEmissionVetoHook.get())
0613       fEmissionVetoHook.reset(new PowhegHooks());
0614 
0615     edm::LogInfo("Pythia8Interface") << "Turning on Emission Veto Hook from pythia8 code";
0616     (fUserHooksVector->hooks).push_back(fEmissionVetoHook);
0617   }
0618 
0619   bool PowhegRes = fMasterGen->settings.flag("POWHEGres:calcScales");
0620   if (PowhegRes) {
0621     edm::LogInfo("Pythia8Interface") << "Turning on resonance scale setting from CMSSW Pythia8Interface";
0622     if (!fPowhegResHook.get())
0623       fPowhegResHook.reset(new PowhegResHook());
0624     (fUserHooksVector->hooks).push_back(fPowhegResHook);
0625   }
0626 
0627   bool PowhegBB4L = fMasterGen->settings.flag("POWHEG:bb4l");
0628   if (PowhegBB4L) {
0629     edm::LogInfo("Pythia8Interface") << "Turning on BB4l hook from CMSSW Pythia8Interface";
0630     if (!fPowhegHooksBB4L.get())
0631       fPowhegHooksBB4L.reset(new PowhegHooksBB4L());
0632     (fUserHooksVector->hooks).push_back(fPowhegHooksBB4L);
0633   }
0634 
0635   bool TopRecoilHook1 = fMasterGen->settings.flag("TopRecoilHook:doTopRecoilIn");
0636   if (TopRecoilHook1) {
0637     edm::LogInfo("Pythia8Interface") << "Turning on RecoilToTop hook from Pythia8Interface";
0638     if (!fTopRecoilHook.get())
0639       fTopRecoilHook.reset(new TopRecoilHook());
0640     (fUserHooksVector->hooks).push_back(fTopRecoilHook);
0641   }
0642 
0643   //adapted from main89.cc in pythia8 examples
0644   bool internalMatching = fMasterGen->settings.flag("JetMatching:merge");
0645   bool internalMerging = !(fMasterGen->settings.word("Merging:Process") == "void");
0646 
0647   if (internalMatching && internalMerging) {
0648     throw edm::Exception(edm::errors::Configuration, "Pythia8Interface")
0649         << " Only one jet matching/merging scheme can be used at a time. \n";
0650   }
0651 
0652   if (internalMatching) {
0653     if (!fJetMatchingPy8InternalHook.get())
0654       fJetMatchingPy8InternalHook.reset(new Pythia8::JetMatchingMadgraph);
0655     (fUserHooksVector->hooks).push_back(fJetMatchingPy8InternalHook);
0656   }
0657 
0658   if (internalMerging) {
0659     int scheme = (fMasterGen->settings.flag("Merging:doUMEPSTree") || fMasterGen->settings.flag("Merging:doUMEPSSubt"))
0660                      ? 1
0661                      : ((fMasterGen->settings.flag("Merging:doUNLOPSTree") ||
0662                          fMasterGen->settings.flag("Merging:doUNLOPSSubt") ||
0663                          fMasterGen->settings.flag("Merging:doUNLOPSLoop") ||
0664                          fMasterGen->settings.flag("Merging:doUNLOPSSubtNLO"))
0665                             ? 2
0666                             : 0);
0667     if (!fMergingHook.get())
0668       fMergingHook.reset(new Pythia8::amcnlo_unitarised_interface(scheme));
0669     (fUserHooksVector->hooks).push_back(fMergingHook);
0670   }
0671 
0672   bool biasedTauDecayer = fMasterGen->settings.flag("BiasedTauDecayer:filter");
0673   if (biasedTauDecayer) {
0674     if (!fBiasedTauDecayer.get()) {
0675       Pythia8::Info localInfo = fMasterGen->info;
0676       fBiasedTauDecayer.reset(new BiasedTauDecayer(&localInfo, &(fMasterGen->settings)));
0677     }
0678     std::vector<int> handledParticles;
0679     handledParticles.push_back(15);
0680     fMasterGen->setDecayPtr(fBiasedTauDecayer, handledParticles);
0681   }
0682 
0683   bool resonanceDecayFilter = fMasterGen->settings.flag("ResonanceDecayFilter:filter");
0684   if (resonanceDecayFilter) {
0685     fResonanceDecayFilterHook.reset(new ResonanceDecayFilterHook);
0686     (fUserHooksVector->hooks).push_back(fResonanceDecayFilterHook);
0687   }
0688 
0689   bool PTFilter = fMasterGen->settings.flag("PTFilter:filter");
0690   if (PTFilter) {
0691     fPTFilterHook.reset(new PTFilterHook);
0692     (fUserHooksVector->hooks).push_back(fPTFilterHook);
0693   }
0694 
0695   if (!(fUserHooksVector->hooks).empty() && !UserHooksSet) {
0696     for (auto &fUserHook : fUserHooksVector->hooks) {
0697       fMasterGen->addUserHooksPtr(fUserHook);
0698     }
0699     UserHooksSet = true;
0700   }
0701 
0702   if (!LHEInputFileName.empty()) {
0703     edm::LogInfo("Pythia8Interface") << "Initialize direct pythia8 reading from LHE file " << LHEInputFileName;
0704     edm::LogInfo("Pythia8Interface") << "Some LHE information can be not stored";
0705     fMasterGen->settings.mode("Beams:frameType", 4);
0706     fMasterGen->settings.word("Beams:LHEF", LHEInputFileName);
0707     status = fMasterGen->init();
0708 
0709   } else {
0710     lhaUP.reset(new LHAupLesHouches());
0711     lhaUP->setScalesFromLHEF(fMasterGen->settings.flag("Beams:setProductionScalesFromLHEF"));
0712     lhaUP->loadRunInfo(lheRunInfo());
0713 
0714     if (fJetMatchingHook.get()) {
0715       fJetMatchingHook->init(lheRunInfo());
0716     }
0717 
0718     fMasterGen->settings.mode("Beams:frameType", 5);
0719     fMasterGen->setLHAupPtr(lhaUP);
0720     edm::LogInfo("Pythia8Interface") << "Initializing MasterGen";
0721     status = fMasterGen->init();
0722   }
0723 
0724   //clean up temp file
0725   if (!slhafile_.empty()) {
0726     std::remove(slhafile_.c_str());
0727   }
0728 
0729   if (pythiaPylistVerbosity > 10) {
0730     if (pythiaPylistVerbosity == 11 || pythiaPylistVerbosity == 13)
0731       fMasterGen->settings.listAll();
0732     if (pythiaPylistVerbosity == 12 || pythiaPylistVerbosity == 13)
0733       fMasterGen->particleData.listAll();
0734   }
0735 
0736   // init decayer
0737   fDecayer->settings.flag("ProcessLevel:all", false);  // trick
0738   fDecayer->settings.flag("ProcessLevel:resonanceDecays", true);
0739   edm::LogInfo("Pythia8Interface") << "Initializing Decayer";
0740   status1 = fDecayer->init();
0741 
0742   if (useEvtGen) {
0743     edm::LogInfo("Pythia8Hadronizer") << "Creating and initializing pythia8 EvtGen plugin";
0744     if (!evtgenDecays.get()) {
0745       evtgenDecays.reset(new EvtGenDecays(fMasterGen.get(), evtgenDecFile, evtgenPdlFile));
0746       for (unsigned int i = 0; i < evtgenUserFiles.size(); i++)
0747         evtgenDecays->readDecayFile(evtgenUserFiles.at(i));
0748     }
0749   }
0750 
0751   return (status && status1);
0752 }
0753 
0754 void Pythia8Hadronizer::statistics() {
0755   fMasterGen->stat();
0756 
0757   if (fEmissionVetoHook.get()) {
0758     edm::LogPrint("Pythia8Interface") << "\n"
0759                                       << "Number of ISR vetoed = " << nISRveto;
0760     edm::LogPrint("Pythia8Interface") << "Number of FSR vetoed = " << nFSRveto;
0761   }
0762 
0763   double xsec = fMasterGen->info.sigmaGen();  // cross section in mb
0764   xsec *= 1.0e9;                              // translate to pb (CMS/Gen "convention" as of May 2009)
0765   double err = fMasterGen->info.sigmaErr();   // cross section err in mb
0766   err *= 1.0e9;                               // translate to pb (CMS/Gen "convention" as of May 2009)
0767   runInfo().setInternalXSec(GenRunInfoProduct::XSec(xsec, err));
0768 }
0769 
0770 bool Pythia8Hadronizer::generatePartonsAndHadronize() {
0771   DJR.resize(0);
0772   nME = -1;
0773   nMEFiltered = -1;
0774 
0775   if (fJetMatchingHook.get()) {
0776     fJetMatchingHook->resetMatchingStatus();
0777     fJetMatchingHook->beforeHadronization(lheEvent());
0778   }
0779 
0780   if (!fMasterGen->next())
0781     return false;
0782 
0783   double mergeweight = fMasterGen.get()->info.mergingWeightNLO();
0784   if (fMergingHook.get()) {
0785     mergeweight *= fMergingHook->getNormFactor();
0786   }
0787 
0788   //protect against 0-weight from ckkw or similar
0789   if (std::abs(mergeweight) == 0.) {
0790     event().reset();
0791     return false;
0792   }
0793 
0794   if (fJetMatchingPy8InternalHook.get()) {
0795     const std::vector<double> djrmatch = fJetMatchingPy8InternalHook->getDJR();
0796     //cap size of djr vector to save storage space (keep only up to first 6 elements)
0797     unsigned int ndjr = std::min(djrmatch.size(), std::vector<double>::size_type(6));
0798     for (unsigned int idjr = 0; idjr < ndjr; ++idjr) {
0799       DJR.push_back(djrmatch[idjr]);
0800     }
0801 
0802     nME = fJetMatchingPy8InternalHook->nMEpartons().first;
0803     nMEFiltered = fJetMatchingPy8InternalHook->nMEpartons().second;
0804   }
0805 
0806   if (evtgenDecays.get())
0807     evtgenDecays->decay();
0808 
0809   event() = std::make_unique<HepMC::GenEvent>();
0810   bool py8hepmc = toHepMC.fill_next_event(*(fMasterGen.get()), event().get());
0811 
0812   if (!py8hepmc) {
0813     return false;
0814   }
0815 
0816   // 0th weight is not filled anymore since v8.30x, pushback manually
0817   event()->weights().push_back(fMasterGen->info.weight());
0818 
0819   //add ckkw/umeps/unlops merging weight
0820   if (mergeweight != 1.) {
0821     event()->weights()[0] *= mergeweight;
0822   }
0823 
0824   if (fEmissionVetoHook.get()) {
0825     nISRveto += fEmissionVetoHook->getNISRveto();
0826     nFSRveto += fEmissionVetoHook->getNFSRveto();
0827   }
0828 
0829   //fill additional weights for systematic uncertainties
0830   if (fMasterGen->info.getWeightsDetailedSize() > 0) {
0831     for (const string &key : fMasterGen->info.initrwgt->weightsKeys) {
0832       double wgt = (*fMasterGen->info.weights_detailed)[key];
0833       event()->weights().push_back(wgt);
0834     }
0835   } else if (fMasterGen->info.getWeightsCompressedSize() > 0) {
0836     for (unsigned int i = 0; i < fMasterGen->info.getWeightsCompressedSize(); i++) {
0837       double wgt = fMasterGen->info.getWeightsCompressedValue(i);
0838       event()->weights().push_back(wgt);
0839     }
0840   }
0841 
0842   // fill shower weights
0843   // http://home.thep.lu.se/~torbjorn/pythia82html/Variations.html
0844   if (fMasterGen->info.nWeights() > 1) {
0845     for (int i = 0; i < fMasterGen->info.nWeights(); ++i) {
0846       double wgt = fMasterGen->info.weight(i);
0847       event()->weights().push_back(wgt);
0848     }
0849   }
0850 
0851 #if 0
0852   // VINCIA shower weights
0853   // http://vincia.hepforge.org/current/share/Vincia/htmldoc/VinciaUncertainties.html
0854   if (fvincia.get()) {
0855     event()->weights()[0] *= fvincia->weight(0);
0856     for (int iVar = 1; iVar < fvincia->nWeights(); iVar++) {
0857       event()->weights().push_back(fvincia->weight(iVar));
0858     }
0859   }
0860 
0861   // Retrieve Dire shower weights
0862   if (fDire.get()) {
0863     fDire->weightsPtr->calcWeight(0.);
0864     fDire->weightsPtr->reset();
0865 
0866     //Make sure the base weight comes first
0867     event()->weights()[0] *= fDire->weightsPtr->getShowerWeight("base");
0868 
0869     unordered_map<string, double>::iterator it;
0870     for (it = fDire->weightsPtr->getShowerWeights()->begin(); it != fDire->weightsPtr->getShowerWeights()->end();
0871          it++) {
0872       if (it->first == "base")
0873         continue;
0874       event()->weights().push_back(it->second);
0875     }
0876   }
0877 #endif
0878 
0879   return true;
0880 }
0881 
0882 bool Pythia8Hadronizer::hadronize() {
0883   DJR.resize(0);
0884   nME = -1;
0885   nMEFiltered = -1;
0886   if (LHEInputFileName.empty())
0887     lhaUP->loadEvent(lheEvent());
0888 
0889   if (fJetMatchingHook.get()) {
0890     fJetMatchingHook->resetMatchingStatus();
0891     fJetMatchingHook->beforeHadronization(lheEvent());
0892   }
0893 
0894   bool py8next = fMasterGen->next();
0895 
0896   double mergeweight = fMasterGen.get()->info.mergingWeightNLO();
0897   if (fMergingHook.get()) {
0898     mergeweight *= fMergingHook->getNormFactor();
0899   }
0900 
0901   //protect against 0-weight from ckkw or similar
0902   if (!py8next || std::abs(mergeweight) == 0.) {
0903     lheEvent()->count(lhef::LHERunInfo::kSelected, 1.0, mergeweight);
0904     event().reset();
0905     return false;
0906   }
0907 
0908   if (fJetMatchingPy8InternalHook.get()) {
0909     const std::vector<double> djrmatch = fJetMatchingPy8InternalHook->getDJR();
0910     //cap size of djr vector to save storage space (keep only up to first 6 elements)
0911     unsigned int ndjr = std::min(djrmatch.size(), std::vector<double>::size_type(6));
0912     for (unsigned int idjr = 0; idjr < ndjr; ++idjr) {
0913       DJR.push_back(djrmatch[idjr]);
0914     }
0915 
0916     nME = fJetMatchingPy8InternalHook->nMEpartons().first;
0917     nMEFiltered = fJetMatchingPy8InternalHook->nMEpartons().second;
0918   }
0919 
0920   // update LHE matching statistics
0921   //
0922   lheEvent()->count(lhef::LHERunInfo::kAccepted, 1.0, mergeweight);
0923 
0924   if (evtgenDecays.get())
0925     evtgenDecays->decay();
0926 
0927   event() = std::make_unique<HepMC::GenEvent>();
0928 
0929   bool py8hepmc = toHepMC.fill_next_event(*(fMasterGen.get()), event().get());
0930 
0931   if (!py8hepmc) {
0932     return false;
0933   }
0934 
0935   // 0th weight is not filled anymore since v8.30x, pushback manually
0936   event()->weights().push_back(fMasterGen->info.weight());
0937 
0938   //add ckkw/umeps/unlops merging weight
0939   if (mergeweight != 1.) {
0940     event()->weights()[0] *= mergeweight;
0941   }
0942 
0943   if (fEmissionVetoHook.get()) {
0944     nISRveto += fEmissionVetoHook->getNISRveto();
0945     nFSRveto += fEmissionVetoHook->getNFSRveto();
0946   }
0947 
0948   // fill shower weights
0949   // http://home.thep.lu.se/~torbjorn/pythia82html/Variations.html
0950   if (fMasterGen->info.nWeights() > 1) {
0951     for (int i = 0; i < fMasterGen->info.nWeights(); ++i) {
0952       double wgt = fMasterGen->info.weight(i);
0953       event()->weights().push_back(wgt);
0954     }
0955   }
0956 
0957   return true;
0958 }
0959 
0960 bool Pythia8Hadronizer::residualDecay() {
0961   Event *pythiaEvent = &(fMasterGen->event);
0962 
0963   int NPartsBeforeDecays = pythiaEvent->size();
0964   int NPartsAfterDecays = event().get()->particles_size();
0965 
0966   if (NPartsAfterDecays == NPartsBeforeDecays)
0967     return true;
0968 
0969   bool result = true;
0970 
0971   for (int ipart = NPartsAfterDecays; ipart > NPartsBeforeDecays; ipart--) {
0972     HepMC::GenParticle *part = event().get()->barcode_to_particle(ipart);
0973 
0974     if (part->status() == 1 && (fDecayer->particleData).canDecay(part->pdg_id())) {
0975       fDecayer->event.reset();
0976       Particle py8part(part->pdg_id(),
0977                        93,
0978                        0,
0979                        0,
0980                        0,
0981                        0,
0982                        0,
0983                        0,
0984                        part->momentum().x(),
0985                        part->momentum().y(),
0986                        part->momentum().z(),
0987                        part->momentum().t(),
0988                        part->generated_mass());
0989       HepMC::GenVertex *ProdVtx = part->production_vertex();
0990       py8part.vProd(ProdVtx->position().x(), ProdVtx->position().y(), ProdVtx->position().z(), ProdVtx->position().t());
0991       py8part.tau((fDecayer->particleData).tau0(part->pdg_id()));
0992       fDecayer->event.append(py8part);
0993       int nentries = fDecayer->event.size();
0994       if (!fDecayer->event[nentries - 1].mayDecay())
0995         continue;
0996       fDecayer->next();
0997       int nentries1 = fDecayer->event.size();
0998       if (nentries1 <= nentries)
0999         continue;  //same number of particles, no decays...
1000 
1001       part->set_status(2);
1002 
1003       result = toHepMC.fill_next_event(*(fDecayer.get()), event().get(), -1, true, part);
1004     }
1005   }
1006 
1007   return result;
1008 }
1009 
1010 void Pythia8Hadronizer::finalizeEvent() {
1011   bool lhe = lheEvent() != nullptr;
1012 
1013   // protection against empty weight container
1014   if ((event()->weights()).empty())
1015     (event()->weights()).push_back(1.);
1016 
1017   // now create the GenEventInfo product from the GenEvent and fill
1018   // the missing pieces
1019   eventInfo() = std::make_unique<GenEventInfoProduct>(event().get());
1020 
1021   // in pythia pthat is used to subdivide samples into different bins
1022   // in LHE mode the binning is done by the external ME generator
1023   // which is likely not pthat, so only filling it for Py6 internal mode
1024   if (!lhe) {
1025     eventInfo()->setBinningValues(std::vector<double>(1, fMasterGen->info.pTHat()));
1026   }
1027 
1028   eventInfo()->setDJR(DJR);
1029   eventInfo()->setNMEPartons(nME);
1030   eventInfo()->setNMEPartonsFiltered(nMEFiltered);
1031 
1032   //******** Verbosity ********
1033 
1034   if (maxEventsToPrint > 0 && (pythiaPylistVerbosity || pythiaHepMCVerbosity || pythiaHepMCVerbosityParticles)) {
1035     maxEventsToPrint--;
1036     if (pythiaPylistVerbosity) {
1037       fMasterGen->info.list();
1038       fMasterGen->event.list();
1039     }
1040 
1041     if (pythiaHepMCVerbosity) {
1042       std::cout << "Event process = " << fMasterGen->info.code() << "\n"
1043                 << "----------------------" << std::endl;
1044       event()->print();
1045     }
1046     if (pythiaHepMCVerbosityParticles) {
1047       std::cout << "Event process = " << fMasterGen->info.code() << "\n"
1048                 << "----------------------" << std::endl;
1049       ascii_io->write_event(event().get());
1050     }
1051   }
1052 }
1053 
1054 std::unique_ptr<GenLumiInfoHeader> Pythia8Hadronizer::getGenLumiInfoHeader() const {
1055   auto genLumiInfoHeader = BaseHadronizer::getGenLumiInfoHeader();
1056 
1057   //fill lhe headers
1058   //*FIXME* initrwgt header is corrupt due to pythia bug
1059   for (const std::string &key : fMasterGen->info.headerKeys()) {
1060     genLumiInfoHeader->lheHeaders().emplace_back(key, fMasterGen->info.header(key));
1061   }
1062 
1063   //check, if it is not only nominal weight
1064   int weights_number = fMasterGen->info.nWeights();
1065   if (fMasterGen->info.initrwgt)
1066     weights_number += fMasterGen->info.initrwgt->weightsKeys.size();
1067   if (weights_number > 1) {
1068     genLumiInfoHeader->weightNames().reserve(weights_number + 1);
1069     genLumiInfoHeader->weightNames().push_back("nominal");
1070   }
1071 
1072   //fill weight names
1073   if (fMasterGen->info.initrwgt) {
1074     for (const std::string &key : fMasterGen->info.initrwgt->weightsKeys) {
1075       std::string weightgroupname;
1076       for (const auto &wgtgrp : fMasterGen->info.initrwgt->weightgroups) {
1077         const auto &wgtgrpwgt = wgtgrp.second.weights.find(key);
1078         if (wgtgrpwgt != wgtgrp.second.weights.end()) {
1079           weightgroupname = wgtgrp.first;
1080         }
1081       }
1082 
1083       std::ostringstream weightname;
1084       weightname << "LHE, id = " << key << ", ";
1085       if (!weightgroupname.empty()) {
1086         weightname << "group = " << weightgroupname << ", ";
1087       }
1088       weightname << fMasterGen->info.initrwgt->weights[key].contents;
1089       genLumiInfoHeader->weightNames().push_back(weightname.str());
1090     }
1091   }
1092 
1093   //fill shower labels
1094   // http://home.thep.lu.se/~torbjorn/pythia82html/Variations.html
1095   // http://home.thep.lu.se/~torbjorn/doxygen/classPythia8_1_1Info.html
1096   if (fMasterGen->info.nWeights() > 1) {
1097     for (int i = 0; i < fMasterGen->info.nWeights(); ++i) {
1098       genLumiInfoHeader->weightNames().push_back(fMasterGen->info.weightLabel(i));
1099     }
1100   }
1101 
1102 #if 0
1103   // VINCIA shower weights
1104   // http://vincia.hepforge.org/current/share/Vincia/htmldoc/VinciaUncertainties.html
1105   if (fvincia.get()) {
1106     for (int iVar = 0; iVar < fvincia->nWeights(); iVar++) {
1107       genLumiInfoHeader->weightNames().push_back(fvincia->weightLabel(iVar));
1108     }
1109   }
1110 
1111   if (fDire.get()) {
1112     //Make sure the base weight comes first
1113     genLumiInfoHeader->weightNames().push_back("base");
1114 
1115     unordered_map<string, double>::iterator it;
1116     for (it = fDire->weightsPtr->getShowerWeights()->begin(); it != fDire->weightsPtr->getShowerWeights()->end();
1117          it++) {
1118       if (it->first == "base")
1119         continue;
1120       genLumiInfoHeader->weightNames().push_back(it->first);
1121     }
1122   }
1123 #endif
1124 
1125   return genLumiInfoHeader;
1126 }
1127 
1128 typedef edm::GeneratorFilter<Pythia8Hadronizer, ExternalDecayDriver> Pythia8GeneratorFilter;
1129 DEFINE_FWK_MODULE(Pythia8GeneratorFilter);
1130 
1131 typedef edm::HadronizerFilter<Pythia8Hadronizer, ExternalDecayDriver> Pythia8HadronizerFilter;
1132 DEFINE_FWK_MODULE(Pythia8HadronizerFilter);
1133 
1134 typedef edm::ConcurrentGeneratorFilter<Pythia8Hadronizer, ConcurrentExternalDecayDriver>
1135     Pythia8ConcurrentGeneratorFilter;
1136 DEFINE_FWK_MODULE(Pythia8ConcurrentGeneratorFilter);
1137 
1138 typedef edm::ConcurrentHadronizerFilter<Pythia8Hadronizer, ConcurrentExternalDecayDriver>
1139     Pythia8ConcurrentHadronizerFilter;
1140 DEFINE_FWK_MODULE(Pythia8ConcurrentHadronizerFilter);