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