Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-07-19 00:58:39

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