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