Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2025-07-04 23:43:57 

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