Pythia8Interface/src/Py8InterfaceBase.cc

0001 #include <memory>
0002 #include <string>
0003
0004 #include "GeneratorInterface/Pythia8Interface/interface/Py8InterfaceBase.h"
0005 #include "GeneratorInterface/Pythia8Interface/interface/SLHAReaderBase.h"
0006
0007 #include "FWCore/Utilities/interface/Exception.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 #include <filesystem>
0010
0011 // EvtGen plugin
0012 //
0013 //#include "Pythia8Plugins/EvtGen.h"
0014
0015 using namespace Pythia8;
0016
0017 namespace gen {
0018
0019   Py8InterfaceBase::Py8InterfaceBase(edm::ParameterSet const& ps)
0020       : BaseHadronizer(ps), useEvtGen(false), evtgenDecays(nullptr) {
0021     fParameters = ps;
0022
0023     pythiaPylistVerbosity = ps.getUntrackedParameter<int>("pythiaPylistVerbosity", 0);
0024     pythiaHepMCVerbosity = ps.getUntrackedParameter<bool>("pythiaHepMCVerbosity", false);
0025     pythiaHepMCVerbosityParticles = ps.getUntrackedParameter<bool>("pythiaHepMCVerbosityParticles", false);
0026     maxEventsToPrint = ps.getUntrackedParameter<int>("maxEventsToPrint", 0);
0027     p8RndmEngine_ = std::make_shared<P8RndmEngine>();
0028
0029     if (pythiaHepMCVerbosityParticles)
0030       ascii_io = new HepMC::IO_AsciiParticles("cout", std::ios::out);
0031
0032     if (ps.exists("useEvtGenPlugin")) {
0033       useEvtGen = true;
0034       auto env = std::getenv("EVTGENDATA");
0035       if (not env) {
0036         throw cms::Exception("EvtGenMissingEnv") << "The environment variable EVTGENDATA must be defined";
0037       }
0038       string evtgenpath(env);
0039       evtgenDecFile = evtgenpath + string("/DECAY_2010.DEC");
0040       evtgenPdlFile = evtgenpath + string("/evt.pdl");
0041
0042       if (ps.exists("evtgenDecFile")) {
0043         edm::FileInPath decay_table(ps.getParameter<std::string>("evtgenDecFile"));
0044         evtgenDecFile = decay_table.fullPath();
0045       }
0046
0047       if (ps.exists("evtgenPdlFile")) {
0048         edm::FileInPath pdt(ps.getParameter<std::string>("evtgenPdlFile"));
0049         evtgenPdlFile = pdt.fullPath();
0050       }
0051
0052       if (ps.exists("evtgenUserFile")) {
0053         std::vector<std::string> user_decays = ps.getParameter<std::vector<std::string> >("evtgenUserFile");
0054         for (unsigned int i = 0; i < user_decays.size(); i++) {
0055           edm::FileInPath user_decay(user_decays.at(i));
0056           evtgenUserFiles.push_back(user_decay.fullPath());
0057         }
0058         //evtgenUserFiles = ps.getParameter< std::vector<std::string> >("evtgenUserFile");
0059       }
0060
0061       if (ps.exists("evtgenUserFileEmbedded")) {
0062         std::vector<std::string> user_decay_lines =
0063             ps.getParameter<std::vector<std::string> >("evtgenUserFileEmbedded");
0064         char tempslhaname[] = "pythia8evtgenXXXXXX";
0065         int fd = mkstemp(tempslhaname);
0066
0067         for (unsigned int i = 0; i < user_decay_lines.size(); i++) {
0068           user_decay_lines.at(i) += "\n";
0069           write(fd, user_decay_lines.at(i).c_str(), user_decay_lines.at(i).size());
0070         }
0071         close(fd);
0072         evtgenUserFiles.push_back(std::string(tempslhaname));
0073       }
0074     }
0075   }
0076
0077   bool Py8InterfaceBase::readSettings(int) {
0078     //Pythia 8's default value for first argument to constructor
0079     const string xmlDir = "../share/Pythia8/xmldoc";
0080     bool printBanner = true;
0081     if (fParameters.exists("printBanner")) {
0082       printBanner = fParameters.getUntrackedParameter<bool>("printBanner");
0083     }
0084     if (!fMasterGen.get())
0085       fMasterGen = std::make_unique<Pythia>(xmlDir, printBanner);
0086     fDecayer = std::make_unique<Pythia>(xmlDir, printBanner);
0087
0088     //add settings for resonance decay filter
0089     fMasterGen->settings.addFlag("BiasedTauDecayer:filter", false);
0090     fMasterGen->settings.addFlag("BiasedTauDecayer:eDecays", true);
0091     fMasterGen->settings.addFlag("BiasedTauDecayer:muDecays", true);
0092
0093     //add settings for resonance decay filter
0094     fMasterGen->settings.addFlag("ResonanceDecayFilter:filter", false);
0095     fMasterGen->settings.addFlag("ResonanceDecayFilter:exclusive", false);
0096     fMasterGen->settings.addFlag("ResonanceDecayFilter:eMuAsEquivalent", false);
0097     fMasterGen->settings.addFlag("ResonanceDecayFilter:eMuTauAsEquivalent", false);
0098     fMasterGen->settings.addFlag("ResonanceDecayFilter:allNuAsEquivalent", false);
0099     fMasterGen->settings.addFlag("ResonanceDecayFilter:udscAsEquivalent", false);
0100     fMasterGen->settings.addFlag("ResonanceDecayFilter:udscbAsEquivalent", false);
0101     fMasterGen->settings.addFlag("ResonanceDecayFilter:wzAsEquivalent", false);
0102     fMasterGen->settings.addMVec("ResonanceDecayFilter:mothers", std::vector<int>(), false, false, 0, 0);
0103     fMasterGen->settings.addMVec("ResonanceDecayFilter:daughters", std::vector<int>(), false, false, 0, 0);
0104
0105     //add settings for PT filter
0106     fMasterGen->settings.addFlag("PTFilter:filter", false);
0107     fMasterGen->settings.addMode("PTFilter:quarkToFilter", 5, true, true, 3, 6);
0108     fMasterGen->settings.addParm("PTFilter:scaleToFilter", 0.4, true, true, 0.0, 10.);
0109     fMasterGen->settings.addParm("PTFilter:quarkRapidity", 10.0, true, true, 0.0, 10.);
0110     fMasterGen->settings.addParm("PTFilter:quarkPt", -.1, true, true, -.1, 100.);
0111
0112     //add settings for RecoilToTop tool
0113     fMasterGen->settings.addFlag("TopRecoilHook:doTopRecoilIn", false);
0114     fMasterGen->settings.addFlag("TopRecoilHook:useOldDipoleIn", false);
0115     fMasterGen->settings.addFlag("TopRecoilHook:doListIn", false);
0116
0117     //add settings for powheg resonance scale calculation
0118     fMasterGen->settings.addFlag("POWHEGres:calcScales", false);
0119     fMasterGen->settings.addFlag("POWHEG:bb4l", false);
0120     fMasterGen->settings.addFlag("POWHEG:bb4l:FSREmission:veto", false);
0121     fMasterGen->settings.addFlag("POWHEG:bb4l:FSREmission:vetoQED", false);
0122     fMasterGen->settings.addFlag("POWHEG:bb4l:DEBUG", false);
0123     fMasterGen->settings.addFlag("POWHEG:bb4l:ScaleResonance:veto", false);
0124     fMasterGen->settings.addFlag("POWHEG:bb4l:FSREmission:vetoDipoleFrame", false);
0125     fMasterGen->settings.addFlag("POWHEG:bb4l:FSREmission:pTpythiaVeto", false);
0126     fMasterGen->settings.addParm("POWHEG:bb4l:pTminVeto", 10.0, true, true, 0.0, 10.);
0127     fMasterGen->settings.addFlag("POWHEG:bb4l:vetoAllRadtypes", false);
0128
0129     fMasterGen->setRndmEnginePtr(p8RndmEngine_);
0130     fDecayer->setRndmEnginePtr(p8RndmEngine_);
0131
0132     fMasterGen->readString("Next:numberShowEvent = 0");
0133     fDecayer->readString("Next:numberShowEvent = 0");
0134
0135     edm::ParameterSet currentParameters;
0136     if (randomIndex() >= 0) {
0137       std::vector<edm::ParameterSet> randomizedParameters =
0138           fParameters.getParameter<std::vector<edm::ParameterSet> >("RandomizedParameters");
0139       currentParameters = randomizedParameters[randomIndex()];
0140     } else {
0141       currentParameters = fParameters;
0142     }
0143
0144     ParameterCollector pCollector = currentParameters.getParameter<edm::ParameterSet>("PythiaParameters");
0145
0146     for (ParameterCollector::const_iterator line = pCollector.begin(); line != pCollector.end(); ++line) {
0147       if (line->find("Random:") != std::string::npos)
0148         throw cms::Exception("PythiaError") << "Attempted to set random number "
0149                                                "using Pythia commands. Please use "
0150                                                "the RandomNumberGeneratorService."
0151                                             << std::endl;
0152
0153       if (!fMasterGen->readString(*line))
0154         throw cms::Exception("PythiaError") << "Pythia 8 did not accept \"" << *line << "\"." << std::endl;
0155
0156       if (line->find("ParticleDecays:") != std::string::npos) {
0157         if (!fDecayer->readString(*line))
0158           throw cms::Exception("PythiaError") << "Pythia 8 Decayer did not accept \"" << *line << "\"." << std::endl;
0159       }
0160     }
0161
0162     slhafile_.clear();
0163
0164     if (currentParameters.exists("SLHAFileForPythia8")) {
0165       std::string slhafilenameshort = currentParameters.getParameter<std::string>("SLHAFileForPythia8");
0166       edm::FileInPath f1(slhafilenameshort);
0167
0168       fMasterGen->settings.mode("SLHA:readFrom", 2);
0169       fMasterGen->settings.word("SLHA:file", f1.fullPath());
0170     } else if (currentParameters.exists("SLHATableForPythia8")) {
0171       std::string slhatable = currentParameters.getParameter<std::string>("SLHATableForPythia8");
0172
0173       makeTmpSLHA(slhatable);
0174     } else if (currentParameters.exists("SLHATreeForPythia8")) {
0175       auto slhaReaderParams = currentParameters.getParameter<edm::ParameterSet>("SLHATreeForPythia8");
0176       std::unique_ptr<SLHAReaderBase> reader =
0177           SLHAReaderFactory::get()->create(slhaReaderParams.getParameter<std::string>("name"), slhaReaderParams);
0178
0179       makeTmpSLHA(reader->getSLHA(currentParameters.getParameter<std::string>("ConfigDescription")));
0180     }
0181
0182     return true;
0183   }
0184
0185   void Py8InterfaceBase::makeTmpSLHA(const std::string& slhatable) {
0186     char tempslhaname[] = "pythia8SLHAtableXXXXXX";
0187     int fd = mkstemp(tempslhaname);
0188     write(fd, slhatable.c_str(), slhatable.size());
0189     close(fd);
0190
0191     slhafile_ = tempslhaname;
0192
0193     fMasterGen->settings.mode("SLHA:readFrom", 2);
0194     fMasterGen->settings.word("SLHA:file", slhafile_);
0195   }
0196
0197   bool Py8InterfaceBase::declareStableParticles(const std::vector<int>& pdgIds) {
0198     for (size_t i = 0; i < pdgIds.size(); i++) {
0199       // FIXME: need to double check if PID's are the same in Py6 & Py8,
0200       //        because the HepPDT translation tool is actually for **Py6**
0201       //
0202       // well, actually it looks like Py8 operates in PDT id's rather than Py6's
0203       //
0204       //    int PyID = HepPID::translatePDTtoPythia( pdgIds[i] );
0205       int PyID = pdgIds[i];
0206       std::ostringstream pyCard;
0207       pyCard << PyID << ":mayDecay=false";
0208
0209       if (fMasterGen->particleData.isParticle(PyID)) {
0210         fMasterGen->readString(pyCard.str());
0211       } else {
0212         edm::LogWarning("DataNotUnderstood") << "Pythia8 does not "
0213                                              << "recognize particle id = " << PyID << std::endl;
0214       }
0215       // alternative:
0216       // set the 2nd input argument warn=false
0217       // - this way Py8 will NOT print warnings about unknown particle code(s)
0218       // fMasterPtr->readString( pyCard.str(), false )
0219     }
0220
0221     return true;
0222   }
0223
0224   bool Py8InterfaceBase::declareSpecialSettings(const std::vector<std::string>& settings) {
0225     for (unsigned int iss = 0; iss < settings.size(); iss++) {
0226       if (settings[iss].find("QED-brem-off") != std::string::npos) {
0227         fMasterGen->readString("TimeShower:QEDshowerByL=off");
0228       } else {
0229         size_t fnd1 = settings[iss].find("Pythia8:");
0230         if (fnd1 != std::string::npos) {
0231           std::string value = settings[iss].substr(fnd1 + 8);
0232           fDecayer->readString(value);
0233         }
0234       }
0235     }
0236     return true;
0237   }
0238
0239   void Py8InterfaceBase::statistics() {
0240     fMasterGen->stat();
0241     return;
0242   }
0243
0244 }  // namespace gen