CMSSW Reference Manual

00001 
00008 #include <iostream>
00009 #include <sstream>
00010 #include <cstdlib>
00011 #include <memory>
00012 #include <cmath>
00013 #include <stdlib.h>
00014 
00015 #include <boost/filesystem.hpp>
00016 
00017 #include <HepMC/GenEvent.h>
00018 #include <HepMC/PdfInfo.h>
00019 #include <HepMC/IO_ExtendedAscii.h>
00020 
00021 #include <ThePEG/Utilities/DynamicLoader.h>
00022 #include <ThePEG/Repository/Repository.h>
00023 #include <ThePEG/Handlers/EventHandler.h>
00024 #include <ThePEG/Handlers/XComb.h>
00025 #include <ThePEG/EventRecord/Event.h>
00026 #include <ThePEG/EventRecord/Particle.h> 
00027 #include <ThePEG/EventRecord/Collision.h>
00028 #include <ThePEG/Config/ThePEG.h>
00029 #include <ThePEG/PDF/PartonExtractor.h>
00030 #include <ThePEG/PDF/PDFBase.h>
00031 
00032 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00033 
00034 #include "GeneratorInterface/ThePEGInterface/interface/ThePEGInterface.h"
00035 #include "GeneratorInterface/ThePEGInterface/interface/HepMCConverter.h"
00036 
00037 using namespace std;
00038 
00039 ThePEGInterface::ThePEGInterface(const edm::ParameterSet &pset) :
00040         dataLocation_(resolveEnvVars(pset.getParameter<string>("dataLocation"))),
00041         generator_(pset.getParameter<string>("generatorModule")),
00042         run_(pset.getParameter<string>("run")),
00043         skipEvents_(pset.getUntrackedParameter<int>("skipEvents", 0))
00044 {
00045         // Write events in hepmc ascii format for debugging purposes
00046         string eventLog = pset.getUntrackedParameter<string>("printEvents", "");
00047         if (!eventLog.empty()) {
00048                 iobc_.reset(new HepMC::IO_ExtendedAscii(eventLog.c_str(), ios::out));
00049                 edm::LogInfo("ThePEGSource") << "Event logging switched on (=> " << eventLog << ")";
00050         }
00051 }
00052 
00053 ThePEGInterface::~ThePEGInterface()
00054 {
00055         if (eg_)
00056                 eg_->finalize();
00057         edm::LogInfo("ThePEGInterface") << "Event generator finalized";
00058 }
00059 
00060 string ThePEGInterface::dataFile(const string &fileName) const
00061 {
00062         return dataLocation_ + "/" + fileName;
00063 }
00064 
00065 string ThePEGInterface::dataFile(const edm::ParameterSet &pset,
00066                                  const string &paramName) const
00067 {
00068         const edm::Entry &entry = pset.retrieve(paramName);
00069         if (entry.typeCode() == 'F')
00070                 return entry.getFileInPath().fullPath();
00071         else
00072                 return dataFile(entry.getString());
00073 }
00074 
00075 string ThePEGInterface::resolveEnvVars(const string &s)
00076 {
00077         string result(s);
00078 
00079         for(;;) {
00080                 string::size_type pos = result.find("${");
00081                 if (pos == string::npos)
00082                         break;
00083 
00084                 string::size_type endpos = result.find('}', pos);
00085                 if (endpos == string::npos)
00086                         break;
00087                 else
00088                         ++endpos;
00089 
00090                 string var = result.substr(pos + 2, endpos - pos - 3);
00091                 const char *path = getenv(var.c_str());
00092 
00093                 result.replace(pos, endpos - pos, path ? path : "");
00094         }
00095 
00096         return result;
00097 }
00098 
00099 void ThePEGInterface::readParameterSet(const edm::ParameterSet &pset, const string &paramSet) const
00100 {
00101         stringstream logstream;
00102 
00103         // Read CMSSW config file parameter set
00104         vector<string> params = pset.getParameter<vector<string> >(paramSet);
00105 
00106         // Loop over the parameter sets
00107         for(vector<string>::const_iterator psIter = params.begin();
00108             psIter != params.end(); ++psIter) {
00109 
00110                 // Include other parameter sets specified by +psName
00111                 if (psIter->find_first_of('+') == 0) {
00112                         edm::LogInfo("ThePEGInterface") << "Loading parameter set (" << psIter->substr(1) << ")";
00113                         readParameterSet(pset, psIter->substr(1));
00114                 }
00115                 // Topmost parameter set is called "parameterSets"
00116                 else if (paramSet == "parameterSets") {
00117                         edm::LogInfo("ThePEGInterface") << "Loading parameter set (" << *psIter << ")";
00118                         readParameterSet(pset, *psIter);
00119                 }
00120                 // Transfer parameters to the repository
00121                 else {
00122                         string line = resolveEnvVars(*psIter);
00123                         string out = ThePEG::Repository::exec(line, logstream);
00124                         if (out != "")
00125                                 edm::LogInfo("ThePEGInterface") << line << " => " << out;
00126                 }
00127         }
00128 }
00129 
00130 void ThePEGInterface::initRepository(const edm::ParameterSet &pset) const
00131 {
00132         /* Initialize the repository from
00133          * 1. the repository file (default values)
00134          * 2. the ThePEG config files
00135          * 3. the CMSSW config blocks
00136          */
00137         stringstream logstream;
00138 
00139         // Read the repository of serialized default values
00140         string repository = dataFile(pset, "repository");
00141         if (!repository.empty()) {
00142                 edm::LogInfo("ThePEGInterface") << "Loading repository (" << repository << ")";
00143                 ThePEG::Repository::load(repository);
00144         }
00145 
00146         if (!getenv("ThePEG_INSTALL_PATH")) {
00147                 vector<string> libdirlist = ThePEG::DynamicLoader::allPaths();
00148                 for(vector<string>::const_iterator libdir = libdirlist.begin();
00149                     libdir < libdirlist.end(); ++libdir) {
00150                         if (libdir->empty() || (*libdir)[0] != '/')
00151                                 continue;
00152                         if (boost::filesystem::exists(*libdir +
00153                                         "/../../share/ThePEG/PDFsets.index")) {
00154                                 setenv("ThePEG_INSTALL_PATH",
00155                                        libdir->c_str(), 0);
00156                                 break;
00157                         }
00158                 }
00159         }
00160 
00161         // Read ThePEG config files to read
00162         vector<string> configFiles = pset.getParameter<vector<string> >("configFiles");
00163 
00164         // Loop over the config files
00165         for(vector<string>::const_iterator iter = configFiles.begin();
00166             iter != configFiles.end(); ++iter) {
00167                 edm::LogInfo("ThePEGInterface") << "Reading config file (" << *iter << ")";
00168                 ThePEG::Repository::read(dataFile(*iter), logstream);
00169                 edm::LogInfo("ThePEGSource") << logstream.str();
00170         }
00171 
00172         // Read CMSSW config file parameter sets starting from "parameterSets"
00173         readParameterSet(pset, "parameterSets");
00174 
00175         // Print the directories where ThePEG looks for libs
00176         vector<string> libdirlist = ThePEG::DynamicLoader::allPaths();
00177         for(vector<string>::const_iterator libdir = libdirlist.begin();
00178             libdir < libdirlist.end(); ++libdir)
00179                 edm::LogInfo("ThePEGInterface") << "DynamicLoader path = " << *libdir << endl;
00180 
00181         // Output status information about the repository
00182         ThePEG::Repository::stats(logstream);
00183         edm::LogInfo("ThePEGInterface") << logstream.str();
00184 }
00185 
00186 void ThePEGInterface::initGenerator()
00187 {
00188         // Get generator from the repository and initialize it
00189         ThePEG::BaseRepository::CheckObjectDirectory(generator_);
00190         ThePEG::EGPtr tmp = ThePEG::BaseRepository::GetObject<ThePEG::EGPtr>(generator_);
00191         if (tmp) {
00192                 eg_ = ThePEG::Repository::makeRun(tmp, run_);
00193                 eg_->initialize();
00194                 edm::LogInfo("ThePEGInterface") << "EventGenerator initialized";
00195         } else
00196                 throw cms::Exception("ThePEGInterface")
00197                         << "EventGenerator could not be initialized!" << endl;
00198 
00199         // Skip events
00200         for (int i = 0; i < skipEvents_; i++) {
00201                 eg_->shoot();
00202                 edm::LogInfo("ThePEGInterface") << "Event discarded";
00203         }
00204 }
00205 
00206 auto_ptr<HepMC::GenEvent> ThePEGInterface::convert(
00207                                         const ThePEG::EventPtr &event)
00208 {
00209         return std::auto_ptr<HepMC::GenEvent>(
00210                 ThePEG::HepMCConverter<HepMC::GenEvent>::convert(*event));
00211 }
00212 
00213 void ThePEGInterface::clearAuxiliary(HepMC::GenEvent *hepmc,
00214                                      HepMC::PdfInfo *pdf)
00215 {
00216         if (hepmc) {
00217                 hepmc->set_event_scale(-1.0);
00218                 hepmc->set_alphaQCD(-1.0);
00219                 hepmc->set_alphaQED(-1.0);
00220         }
00221         if (pdf) {
00222                 pdf->set_id1(-100);
00223                 pdf->set_id2(-100);
00224                 pdf->set_x1(-1.0);
00225                 pdf->set_x2(-1.0);
00226                 pdf->set_scalePDF(-1.0);
00227                 pdf->set_pdf1(-1.0);
00228                 pdf->set_pdf2(-1.0);
00229         }
00230 }
00231 
00232 void ThePEGInterface::fillAuxiliary(HepMC::GenEvent *hepmc,
00233                                     HepMC::PdfInfo *pdf,
00234                                     const ThePEG::EventPtr &event)
00235 {
00236         if (!event->primaryCollision())
00237                 return;
00238 
00239         ThePEG::tcEHPtr eh = ThePEG::dynamic_ptr_cast<ThePEG::tcEHPtr>(
00240                                 event->primaryCollision()->handler());
00241         double scale = eh->lastScale();
00242 
00243         if (hepmc) {
00244                 if (hepmc->event_scale() < 0 && scale > 0)
00245                         hepmc->set_event_scale(std::sqrt(scale) / ThePEG::GeV);
00246 
00247                 if (hepmc->alphaQCD() < 0)
00248                         hepmc->set_alphaQCD(eh->lastAlphaS());
00249                 if (hepmc->alphaQED() < 0)
00250                         hepmc->set_alphaQED(eh->lastAlphaEM());
00251         }
00252 
00253         if (pdf) {
00254                 const ThePEG::PPair &beams = eh->lastParticles();
00255                 const ThePEG::PPair &partons = eh->lastPartons();
00256                 ThePEG::tcPDFPtr pdf1 = eh->lastExtractor()->getPDF(
00257                                                 beams.first->dataPtr());
00258                 ThePEG::tcPDFPtr pdf2 = eh->lastExtractor()->getPDF(
00259                                                 beams.second->dataPtr());
00260                 double x1, x2;
00261 
00262                 if (pdf->id1() == -100) {
00263                         int id = partons.first->id();
00264                         pdf->set_id1(id == 21 ? 0 : id);
00265                 }
00266                 if (pdf->id2() == -100) {
00267                         int id = partons.second->id();
00268                         pdf->set_id2(id == 21 ? 0 : id);
00269                 }
00270 
00271                 if (pdf->scalePDF() < 0)
00272                         pdf->set_scalePDF(std::sqrt(scale) / ThePEG::GeV);
00273                 else
00274                         scale = ThePEG::sqr(pdf->scalePDF()) * ThePEG::GeV;
00275 
00276                 if (pdf->x1() < 0) {
00277                         x1 = eh->lastX1();
00278                         pdf->set_x1(x1);
00279                 } else
00280                         x1 = pdf->x1();
00281 
00282                 if (pdf->x2() < 0) {
00283                         x2 = eh->lastX2();
00284                         pdf->set_x2(x2);
00285                 } else
00286                         x2 = pdf->x2();
00287 
00288                 if (pdf1 && pdf->pdf1() < 0) {
00289                         double v = pdf1->xfx(beams.first->dataPtr(),
00290                                              partons.first->dataPtr(),
00291                                              scale, x1);
00292                         if (v > 0)
00293                                 pdf->set_pdf1(v);
00294                         else
00295                                 pdf->set_pdf2(-1.0);
00296                 }
00297                 if (pdf2 && pdf->pdf2() < 0) {
00298                         double v = pdf2->xfx(beams.first->dataPtr(),
00299                                              partons.first->dataPtr(),
00300                                              scale, x2);
00301                         if (v > 0)
00302                                 pdf->set_pdf2(v);
00303                         else
00304                                 pdf->set_pdf2(-1.0);
00305                 }
00306         }
00307 
00308 }
ThePEGInterface.cc
