Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-04-04 01:26:40 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!