Pythia6Interface/src/Pythia6Service.cc

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 // #include "GeneratorInterface/Core/interface/ParameterCollector.h"
0024
0025 // This will force the symbols below to be kept, even in the case pythia6
0026 // is an archive library.
0027 //extern "C" void pyexec_(void);
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   // getInstance will throw if no one used enter/leave
0061   // or this is the wrong caller class, like e.g. Herwig6Instance
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    ParameterCollector collector(ps.getParameter<edm::ParameterSet>("PythiaParameters"));
0082
0083    fParamGeneral.clear();
0084    fParamCSA.clear();
0085    fParamSLHA.clear();
0086
0087    fParamGeneral = std::vector<std::string>(collector.begin(), collector.end());
0088    fParamCSA     = std::vector<std::string>(collector.begin("CSAParameters"), collector.end());
0089    fParamSLHA    = std::vector<std::string>(collector.begin("SLHAParameters"), collector.end());
0090 */
0091
0092   // Set PYTHIA parameters in a single ParameterSet
0093   //
0094   edm::ParameterSet pythia_params = ps.getParameter<edm::ParameterSet>("PythiaParameters");
0095
0096   // read and sort Pythia6 cards
0097   //
0098   std::vector<std::string> setNames = pythia_params.getParameter<std::vector<std::string> >("parameterSets");
0099
0100   // std::vector<std::string>   paramLines;
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   // now pass general config cards
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     // Null pad the string should not be needed because it uses
0174     // read, which will look for \n, but just in case...
0175     for (size_t i = 0; i < SETCSAPARBUFSIZE; ++i)
0176       buf[i] = ' ';
0177     // Skip empty parameters.
0178     if (iter->length() <= 0)
0179       continue;
0180     // Limit the size of the string to something which fits the buffer.
0181     size_t maxSize = iter->length() > (SETCSAPARBUFSIZE - 2) ? (SETCSAPARBUFSIZE - 2) : iter->length();
0182     strncpy(buf, iter->c_str(), maxSize);
0183     // Add extra \n if missing, otherwise "read" continues reading.
0184     if (buf[maxSize - 1] != '\n') {
0185       buf[maxSize] = '\n';
0186       // Null terminate in case the string is passed back to C.
0187       // Not sure that is actually needed.
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     // Write Pythia particle table to this card file.
0215     call_pyupda(1, fUnitPYUPDA);
0216   } else {
0217     std::cout << "=== READING PYUPDA FILE " << file << " ===" << std::endl;
0218     fioopn_(&fUnitPYUPDA, file, strlen(file));
0219     // Update Pythia particle table with this card file.
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     //std::cout << " start, end = " << start << " " << end << std::endl;
0252     std::string shortfile = iter->substr(start, end - start + 1);
0253     FileInPath f1(shortfile);
0254     std::string file = f1.fullPath();
0255
0256     /*
0257     //
0258     // override what might have be given via the external config
0259     //
0260         std::ostringstream pyCard ;
0261         pyCard << "IMSS(21)=" << fUnitSLHA;
0262         call_pygive( pyCard.str() );
0263         pyCard << "IMSS(22)=" << fUnitSLHA;
0264         call_pygive( pyCard.str() );
0265
0266     fioopn_( &fUnitSLHA, file.c_str(), file.length() );
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   //   std::cout<<"=== CALLING setPYUPDAParams === "<<afterPyinit<<" "<<fParamPYUPDA.size()<<std::endl;
0281
0282   // This assumes that PYUPDAFILE only appears once ...
0283
0284   for (std::vector<std::string>::const_iterator iter = fParamPYUPDA.begin(); iter != fParamPYUPDA.end(); iter++) {
0285     //     std::cout<<"PYUPDA check "<<*iter<<std::endl;
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       //std::cout << " start, end = " << start << " " << end << std::endl;
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       // If reading, get full path to file and require it to exist.
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     int unit = 24;
0398     fioopn_(&unit, fname, std::strlen(fname));
0399     std::remove(fname);
0400
0401     call_pygive("IMSS(21)=24");
0402     call_pygive("IMSS(22)=24");
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 }