File indexing completed on 2025-01-09 23:33:35
0001 #include <algorithm>
0002 #include <functional>
0003 #include <iostream>
0004 #include <sstream>
0005 #include <fstream>
0006 #include <cmath>
0007 #include <string>
0008 #include <set>
0009
0010 #include <boost/algorithm/string/classification.hpp>
0011 #include <boost/algorithm/string/split.hpp>
0012 #include <filesystem>
0013
0014 #include "CLHEP/Random/RandomEngine.h"
0015
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017
0018 #include "FWCore/ParameterSet/interface/FileInPath.h"
0019 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0020
0021 #include "GeneratorInterface/Pythia6Interface/interface/Pythia6Service.h"
0022 #include "GeneratorInterface/Pythia6Interface/interface/Pythia6Declarations.h"
0023
0024
0025
0026
0027
0028 extern "C" void pyedit_(void);
0029 __attribute__((visibility("hidden"))) void dummy() {
0030 using namespace gen;
0031 int dummy = 0;
0032 double dummy2 = 0;
0033 char* dummy3 = nullptr;
0034 pyexec_();
0035 pystat_(nullptr);
0036 pyjoin_(dummy, &dummy);
0037 py1ent_(dummy, dummy, dummy2, dummy2, dummy2);
0038 pygive_(dummy3, dummy);
0039 pycomp_(dummy);
0040 pylist_(nullptr);
0041 pyevnt_();
0042 pyedit_();
0043 }
0044
0045 extern "C" {
0046 void fioopn_(int* unit, const char* line, int length);
0047 void fioopnw_(int* unit, const char* line, int length);
0048 void fiocls_(int* unit);
0049 void pyslha_(int*, int*, int*);
0050 static int call_pyslha(int mupda, int kforig = 0) {
0051 int iretrn = 0;
0052 pyslha_(&mupda, &kforig, &iretrn);
0053 return iretrn;
0054 }
0055
0056 void pyupda_(int*, int*);
0057 static void call_pyupda(int opt, int iunit) { pyupda_(&opt, &iunit); }
0058
0059 double gen::pyr_(int* idummy) {
0060
0061
0062 Pythia6Service* service = FortranInstance::getInstance<Pythia6Service>();
0063 return service->fRandomEngine->flat();
0064 }
0065 }
0066
0067 using namespace gen;
0068 using namespace edm;
0069
0070 Pythia6Service* Pythia6Service::fPythia6Owner = nullptr;
0071
0072 Pythia6Service::Pythia6Service() : fRandomEngine(nullptr), fUnitSLHA(24), fUnitPYUPDA(25) {}
0073
0074 Pythia6Service::Pythia6Service(const ParameterSet& ps) : fRandomEngine(nullptr), fUnitSLHA(24), fUnitPYUPDA(25) {
0075 if (fPythia6Owner)
0076 throw cms::Exception("PythiaError") << "Two Pythia6Service instances claiming Pythia6 ownership." << std::endl;
0077
0078 fPythia6Owner = this;
0079
0080
0081
0082
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092
0093
0094 edm::ParameterSet pythia_params = ps.getParameter<edm::ParameterSet>("PythiaParameters");
0095
0096
0097
0098 std::vector<std::string> setNames = pythia_params.getParameter<std::vector<std::string> >("parameterSets");
0099
0100
0101 fParamGeneral.clear();
0102 fParamCSA.clear();
0103 fParamSLHA.clear();
0104 fParamPYUPDA.clear();
0105
0106 for (std::vector<std::string>::const_iterator iter = setNames.begin(); iter != setNames.end(); ++iter) {
0107 std::vector<std::string> lines = pythia_params.getParameter<std::vector<std::string> >(*iter);
0108
0109 for (std::vector<std::string>::const_iterator line = lines.begin(); line != lines.end(); ++line) {
0110 if (line->substr(0, 7) == "MRPY(1)")
0111 throw cms::Exception("PythiaError") << "Attempted to set random number"
0112 " using Pythia command 'MRPY(1)'."
0113 " Please use the"
0114 " RandomNumberGeneratorService."
0115 << std::endl;
0116
0117 if (*iter == "CSAParameters") {
0118 fParamCSA.push_back(*line);
0119 } else if (*iter == "SLHAParameters") {
0120 fParamSLHA.push_back(*line);
0121 } else if (*iter == "PYUPDAParameters") {
0122 fParamPYUPDA.push_back(*line);
0123 } else {
0124 fParamGeneral.push_back(*line);
0125 }
0126 }
0127 }
0128 }
0129
0130 Pythia6Service::~Pythia6Service() {
0131 if (fPythia6Owner == this)
0132 fPythia6Owner = nullptr;
0133
0134 fParamGeneral.clear();
0135 fParamCSA.clear();
0136 fParamSLHA.clear();
0137 fParamPYUPDA.clear();
0138 }
0139
0140 void Pythia6Service::enter() {
0141 FortranInstance::enter();
0142
0143 if (!fPythia6Owner) {
0144 edm::LogInfo("Generator|Pythia6Interface") << "gen::Pythia6Service is going to initialise Pythia, as no other "
0145 "instace has done so yet, and Pythia service routines have been "
0146 "requested by a dummy instance."
0147 << std::endl;
0148
0149 call_pygive("MSTU(12)=12345");
0150 call_pyinit("NONE", "", "", 0.0);
0151
0152 fPythia6Owner = this;
0153 }
0154 }
0155
0156 void Pythia6Service::setGeneralParams() {
0157
0158
0159 for (std::vector<std::string>::const_iterator iter = fParamGeneral.begin(); iter != fParamGeneral.end(); ++iter) {
0160 if (!call_pygive(*iter))
0161 throw cms::Exception("PythiaError") << "Pythia did not accept \"" << *iter << "\"." << std::endl;
0162 }
0163
0164 return;
0165 }
0166
0167 void Pythia6Service::setCSAParams() {
0168 #define SETCSAPARBUFSIZE 514
0169 char buf[SETCSAPARBUFSIZE];
0170
0171 txgive_init_();
0172 for (std::vector<std::string>::const_iterator iter = fParamCSA.begin(); iter != fParamCSA.end(); ++iter) {
0173
0174
0175 for (size_t i = 0; i < SETCSAPARBUFSIZE; ++i)
0176 buf[i] = ' ';
0177
0178 if (iter->empty())
0179 continue;
0180
0181 size_t maxSize = iter->length() > (SETCSAPARBUFSIZE - 2) ? (SETCSAPARBUFSIZE - 2) : iter->length();
0182 strncpy(buf, iter->c_str(), maxSize);
0183
0184 if (buf[maxSize - 1] != '\n') {
0185 buf[maxSize] = '\n';
0186
0187
0188 buf[maxSize + 1] = 0;
0189 }
0190 txgive_(buf, iter->length());
0191 }
0192
0193 return;
0194 #undef SETCSAPARBUFSIZE
0195 }
0196
0197 void Pythia6Service::openSLHA(const char* file) {
0198 std::ostringstream pyCard1;
0199 pyCard1 << "IMSS(21)=" << fUnitSLHA;
0200 call_pygive(pyCard1.str());
0201 std::ostringstream pyCard2;
0202 pyCard2 << "IMSS(22)=" << fUnitSLHA;
0203 call_pygive(pyCard2.str());
0204
0205 fioopn_(&fUnitSLHA, file, strlen(file));
0206
0207 return;
0208 }
0209
0210 void Pythia6Service::openPYUPDA(const char* file, bool write_file) {
0211 if (write_file) {
0212 std::cout << "=== WRITING PYUPDA FILE " << file << " ===" << std::endl;
0213 fioopnw_(&fUnitPYUPDA, file, strlen(file));
0214
0215 call_pyupda(1, fUnitPYUPDA);
0216 } else {
0217 std::cout << "=== READING PYUPDA FILE " << file << " ===" << std::endl;
0218 fioopn_(&fUnitPYUPDA, file, strlen(file));
0219
0220 call_pyupda(3, fUnitPYUPDA);
0221 }
0222
0223 return;
0224 }
0225
0226 void Pythia6Service::closeSLHA() {
0227 fiocls_(&fUnitSLHA);
0228
0229 return;
0230 }
0231
0232 void Pythia6Service::closePYUPDA() {
0233 fiocls_(&fUnitPYUPDA);
0234
0235 return;
0236 }
0237
0238 void Pythia6Service::setSLHAParams() {
0239 for (std::vector<std::string>::const_iterator iter = fParamSLHA.begin(); iter != fParamSLHA.end(); iter++) {
0240 if (iter->find("SLHAFILE", 0) == std::string::npos)
0241 continue;
0242 std::string::size_type start = iter->find_first_of("=") + 1;
0243 std::string::size_type end = iter->length() - 1;
0244 std::string::size_type temp = iter->find_first_of("'", start);
0245 if (temp != std::string::npos) {
0246 start = temp + 1;
0247 end = iter->find_last_of("'") - 1;
0248 }
0249 start = iter->find_first_not_of(" ", start);
0250 end = iter->find_last_not_of(" ", end);
0251
0252 std::string shortfile = iter->substr(start, end - start + 1);
0253 FileInPath f1(shortfile);
0254 const std::string& file = f1.fullPath();
0255
0256
0257
0258
0259
0260
0261
0262
0263
0264
0265
0266
0267
0268
0269 openSLHA(file.c_str());
0270 }
0271
0272 return;
0273 }
0274
0275 void Pythia6Service::setPYUPDAParams(bool afterPyinit) {
0276 std::string shortfile;
0277 bool write_file = false;
0278 bool usePostPyinit = false;
0279
0280
0281
0282
0283
0284 for (std::vector<std::string>::const_iterator iter = fParamPYUPDA.begin(); iter != fParamPYUPDA.end(); iter++) {
0285
0286 if (iter->find("PYUPDAFILE", 0) != std::string::npos) {
0287 std::string::size_type start = iter->find_first_of("=") + 1;
0288 std::string::size_type end = iter->length() - 1;
0289 std::string::size_type temp = iter->find_first_of("'", start);
0290 if (temp != std::string::npos) {
0291 start = temp + 1;
0292 end = iter->find_last_of("'") - 1;
0293 }
0294 start = iter->find_first_not_of(" ", start);
0295 end = iter->find_last_not_of(" ", end);
0296
0297 shortfile = iter->substr(start, end - start + 1);
0298 } else if (iter->find("PYUPDAWRITE", 0) != std::string::npos) {
0299 write_file = true;
0300 } else if (iter->find("PYUPDApostPYINIT", 0) != std::string::npos) {
0301 usePostPyinit = true;
0302 }
0303 }
0304
0305 if (!shortfile.empty()) {
0306 std::string file;
0307 if (write_file) {
0308 file = shortfile;
0309 } else {
0310
0311 FileInPath f1(shortfile);
0312 file = f1.fullPath();
0313 }
0314
0315 if (afterPyinit == usePostPyinit || (write_file && afterPyinit)) {
0316 openPYUPDA(file.c_str(), write_file);
0317 }
0318 }
0319
0320 return;
0321 }
0322
0323 void Pythia6Service::setSLHAFromHeader(const std::vector<std::string>& lines) {
0324 std::set<std::string> blocks;
0325 unsigned int model = 0, subModel = 0;
0326 char tempslhaname[] = "pythia6slhaXXXXXX";
0327 int fd = mkstemp(tempslhaname);
0328
0329 std::string block;
0330 std::stringstream f_info;
0331 for (std::vector<std::string>::const_iterator iter = lines.begin(); iter != lines.end(); ++iter) {
0332 f_info << *iter;
0333
0334 std::string line = *iter;
0335 std::transform(line.begin(), line.end(), line.begin(), (int (*)(int))std::toupper);
0336 std::string::size_type pos = line.find('#');
0337 if (pos != std::string::npos)
0338 line.resize(pos);
0339
0340 if (line.empty())
0341 continue;
0342
0343 if (!boost::algorithm::is_space()(line[0])) {
0344 std::vector<std::string> tokens;
0345 boost::split(tokens, line, boost::algorithm::is_space(), boost::token_compress_on);
0346 if (tokens.empty())
0347 continue;
0348 block.clear();
0349 if (tokens.size() < 2)
0350 continue;
0351 if (tokens[0] == "BLOCK") {
0352 block = tokens[1];
0353 blocks.insert(block);
0354 continue;
0355 }
0356
0357 if (tokens[0] == "DECAY") {
0358 block = "DECAY";
0359 blocks.insert(block);
0360 }
0361 } else if (block == "MODSEL") {
0362 std::istringstream ss(line);
0363 ss >> model >> subModel;
0364 } else if (block == "SMINPUTS") {
0365 std::istringstream ss(line);
0366 int index;
0367 double value;
0368 ss >> index >> value;
0369 switch (index) {
0370 case 1:
0371 pydat1_.paru[103 - 1] = 1.0 / value;
0372 break;
0373 case 2:
0374 pydat1_.paru[105 - 1] = value;
0375 break;
0376 case 4:
0377 pydat2_.pmas[0][23 - 1] = value;
0378 break;
0379 case 6:
0380 pydat2_.pmas[0][6 - 1] = value;
0381 break;
0382 case 7:
0383 pydat2_.pmas[0][15 - 1] = value;
0384 break;
0385 }
0386 }
0387 }
0388 write(fd, f_info.str().c_str(), f_info.str().size());
0389 close(fd);
0390
0391 if (blocks.count("SMINPUTS"))
0392 pydat1_.paru[102 - 1] =
0393 0.5 - std::sqrt(0.25 - pydat1_.paru[0] * M_SQRT1_2 * pydat1_.paru[103 - 1] / pydat1_.paru[105 - 1] /
0394 (pydat2_.pmas[0][23 - 1] * pydat2_.pmas[0][23 - 1]));
0395
0396
0397
0398
0399
0400
0401
0402
0403
0404
0405 openSLHA(tempslhaname);
0406 std::remove(tempslhaname);
0407
0408 if (model || blocks.count("HIGMIX") || blocks.count("SBOTMIX") || blocks.count("STOPMIX") ||
0409 blocks.count("STAUMIX") || blocks.count("AMIX") || blocks.count("NMIX") || blocks.count("UMIX") ||
0410 blocks.count("VMIX"))
0411 call_pyslha(1);
0412 if (model || blocks.count("QNUMBERS") || blocks.count("PARTICLE") || blocks.count("MINPAR") ||
0413 blocks.count("EXTPAR") || blocks.count("SMINPUTS") || blocks.count("SMINPUTS"))
0414 call_pyslha(0);
0415 if (blocks.count("MASS"))
0416 call_pyslha(5, 0);
0417 if (blocks.count("DECAY"))
0418 call_pyslha(2);
0419
0420 return;
0421 }