LHEInterface/src/LHEReader.cc

0001 #include <algorithm>
0002 #include <iomanip>
0003 #include <iostream>
0004 #include <memory>
0005
0006 #include <cstdio>
0007 #include <cstring>
0008 #include <fstream>
0009 #include <sstream>
0010 #include <string>
0011 #include <vector>
0012
0013 #include <xercesc/sax2/Attributes.hpp>
0014 #include <xercesc/dom/DOM.hpp>
0015
0016 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0017 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0018 #include "FWCore/Utilities/interface/Exception.h"
0019 #include "FWCore/Utilities/interface/TimeOfDay.h"
0020
0021 #include "GeneratorInterface/LHEInterface/interface/LHEReader.h"
0022 #include "GeneratorInterface/LHEInterface/interface/LHERunInfo.h"
0023 #include "GeneratorInterface/LHEInterface/interface/LHEEvent.h"
0024
0025 #include "Utilities/StorageFactory/interface/IOTypes.h"
0026 #include "Utilities/StorageFactory/interface/Storage.h"
0027 #include "Utilities/StorageFactory/interface/StorageFactory.h"
0028
0029 #include "XMLUtils.h"
0030
0031 XERCES_CPP_NAMESPACE_USE
0032
0033 namespace lhef {
0034
0035   static void logFileAction(char const *msg, std::string const &fileName) {
0036     edm::LogAbsolute("fileAction") << std::setprecision(0) << edm::TimeOfDay() << msg << fileName;
0037     edm::FlushMessageLog();
0038   }
0039
0040   class LHEReader::Source {
0041   public:
0042     Source() {}
0043     virtual ~Source() {}
0044     virtual XMLDocument *createReader(XMLDocument::Handler &handler) = 0;
0045   };
0046
0047   class LHEReader::FileSource : public LHEReader::Source {
0048   public:
0049     FileSource(const std::string &fileURL) {
0050       using namespace edm::storage;
0051       auto storage = StorageFactory::get()->open(fileURL, IOFlags::OpenRead);
0052
0053       if (!storage)
0054         throw cms::Exception("FileOpenError")
0055             << "Could not open LHE file \"" << fileURL << "\" for reading" << std::endl;
0056
0057       fileStream = std::make_unique<StorageWrap>(std::move(storage));
0058     }
0059
0060     ~FileSource() override {}
0061
0062     XMLDocument *createReader(XMLDocument::Handler &handler) override { return new XMLDocument(fileStream, handler); }
0063
0064   private:
0065     std::unique_ptr<StorageWrap> fileStream;
0066   };
0067
0068   class LHEReader::StringSource : public LHEReader::Source {
0069   public:
0070     StringSource(const std::string &inputs) {
0071       if (inputs.empty())
0072         throw cms::Exception("StreamOpenError") << "Empty LHE file string name \"" << std::endl;
0073
0074       std::stringstream *tmpis = new std::stringstream(inputs);
0075       fileStream.reset(tmpis);
0076     }
0077
0078     ~StringSource() override {}
0079
0080     XMLDocument *createReader(XMLDocument::Handler &handler) override { return new XMLDocument(fileStream, handler); }
0081
0082   private:
0083     std::unique_ptr<std::istream> fileStream;
0084   };
0085
0086   class LHEReader::XMLHandler : public XMLDocument::Handler {
0087   public:
0088     typedef std::vector<std::pair<std::string, std::string> > wgt_info;
0089     XMLHandler()
0090         : impl(nullptr),
0091           gotObject(kNone),
0092           mode(kNone),
0093           xmlHeader(nullptr),
0094           xmlEvent(nullptr),
0095           headerOk(false),
0096           npLO(-99),
0097           npNLO(-99) {}
0098     ~XMLHandler() override {
0099       if (xmlHeader)
0100         xmlHeader->release();
0101       if (xmlEvent)
0102         xmlEvent->release();
0103     }
0104
0105     enum Object { kNone = 0, kHeader, kInit, kComment, kEvent };
0106
0107     void reset() {
0108       headerOk = false;
0109       weightsinevent.clear();
0110       gotObject = kNone;
0111       mode = kNone;
0112     }
0113
0114     const wgt_info &weightInfo() const { return weightsinevent; }
0115
0116   protected:
0117     void startElement(const XMLCh *const uri,
0118                       const XMLCh *const localname,
0119                       const XMLCh *const qname,
0120                       const Attributes &attributes) override;
0121
0122     void endElement(const XMLCh *const uri, const XMLCh *const localname, const XMLCh *const qname) override;
0123
0124     void characters(const XMLCh *const chars, const XMLSize_t length) override;
0125     void comment(const XMLCh *const chars, const XMLSize_t length) override;
0126
0127   private:
0128     friend class LHEReader;
0129
0130     bool skipEvent = false;
0131     std::unique_ptr<DOMImplementation> impl;
0132     std::string buffer;
0133     Object gotObject;
0134     Object mode;
0135     DOMDocument *xmlHeader;
0136     DOMDocument *xmlEvent;
0137     std::vector<DOMElement *> xmlNodes, xmlEventNodes;
0138     bool headerOk;
0139     std::vector<LHERunInfo::Header> headers;
0140     wgt_info weightsinevent;
0141     int npLO;
0142     int npNLO;
0143     std::vector<float> scales;
0144     int evtnum = -1;
0145   };
0146
0147   static void attributesToDom(DOMElement *dom, const Attributes &attributes) {
0148     for (unsigned int i = 0; i < attributes.getLength(); i++) {
0149       const XMLCh *name = attributes.getQName(i);
0150       const XMLCh *value = attributes.getValue(i);
0151
0152       dom->setAttribute(name, value);
0153     }
0154   }
0155
0156   static void fillHeader(LHERunInfo::Header &header, const char *data, int len = -1) {
0157     const char *end = len >= 0 ? (data + len) : nullptr;
0158     while (*data && (!end || data < end)) {
0159       std::size_t len = std::strcspn(data, "\r\n");
0160       if (end && data + len > end)
0161         len = end - data;
0162       if (data[len] == '\r' && data[len + 1] == '\n')
0163         len += 2;
0164       else if (data[len])
0165         len++;
0166       header.addLine(std::string(data, len));
0167       data += len;
0168     }
0169   }
0170
0171   void LHEReader::XMLHandler::startElement(const XMLCh *const uri,
0172                                            const XMLCh *const localname,
0173                                            const XMLCh *const qname,
0174                                            const Attributes &attributes) {
0175     std::string name((const char *)XMLSimpleStr(qname));
0176
0177     if (!headerOk) {
0178       if (name != "LesHouchesEvents")
0179         throw cms::Exception("InvalidFormat") << "LHE file has invalid header" << std::endl;
0180       headerOk = true;
0181       return;
0182     }
0183
0184     if (mode == kHeader) {
0185       DOMElement *elem = xmlHeader->createElement(qname);
0186       attributesToDom(elem, attributes);
0187       xmlNodes.back()->appendChild(elem);
0188       xmlNodes.push_back(elem);
0189       return;
0190     } else if (mode == kEvent) {
0191       if (skipEvent) {
0192         return;
0193       }
0194
0195       DOMElement *elem = xmlEvent->createElement(qname);
0196       attributesToDom(elem, attributes);
0197
0198       //TODO this is a hack (even more than the rest of this class)
0199       if (name == "rwgt") {
0200         xmlEventNodes[0]->appendChild(elem);
0201       } else if (name == "wgt") {
0202         xmlEventNodes[1]->appendChild(elem);
0203       } else if (name == "scales") {
0204         for (XMLSize_t iscale = 0; iscale < attributes.getLength(); ++iscale) {
0205           int ipart = 0;
0206           const char *scalename = XMLSimpleStr(attributes.getQName(iscale));
0207           int nmatch = sscanf(scalename, "pt_clust_%d", &ipart);
0208
0209           if (nmatch != 1) {
0210             edm::LogError("Generator|LHEInterface") << "invalid attribute in <scales> tag" << std::endl;
0211           }
0212
0213           float scaleval;
0214           const char *scalevalstr = XMLSimpleStr(attributes.getValue(iscale));
0215           sscanf(scalevalstr, "%e", &scaleval);
0216
0217           scales.push_back(scaleval);
0218         }
0219       } else if (name == "event_num") {
0220         const char *evtnumstr = XMLSimpleStr(attributes.getValue(XMLString::transcode("num")));
0221         sscanf(evtnumstr, "%d", &evtnum);
0222       }
0223       xmlEventNodes.push_back(elem);
0224       return;
0225     } else if (mode == kInit) {
0226       //skip unknown tags in init block as well
0227       return;
0228     } else if (mode != kNone) {
0229       throw cms::Exception("InvalidFormat") << "LHE file has invalid format" << std::endl;
0230     }
0231
0232     if (name == "header") {
0233       if (!impl)
0234         impl.reset(DOMImplementationRegistry::getDOMImplementation(XMLUniStr("Core")));
0235
0236       xmlHeader = impl->createDocument(nullptr, qname, nullptr);
0237       xmlNodes.resize(1);
0238       xmlNodes[0] = xmlHeader->getDocumentElement();
0239       mode = kHeader;
0240     }
0241     if (name == "init") {
0242       mode = kInit;
0243     } else if (name == "event") {
0244       if (!skipEvent) {
0245         if (!impl)
0246           impl.reset(DOMImplementationRegistry::getDOMImplementation(XMLUniStr("Core")));
0247
0248         if (xmlEvent)
0249           xmlEvent->release();
0250         xmlEvent = impl->createDocument(nullptr, qname, nullptr);
0251         weightsinevent.resize(0);
0252         scales.clear();
0253
0254         npLO = -99;
0255         npNLO = -99;
0256         const XMLCh *npLOval = attributes.getValue(XMLString::transcode("npLO"));
0257         if (npLOval) {
0258           const char *npLOs = XMLSimpleStr(npLOval);
0259           sscanf(npLOs, "%d", &npLO);
0260         }
0261         const XMLCh *npNLOval = attributes.getValue(XMLString::transcode("npNLO"));
0262         if (npNLOval) {
0263           const char *npNLOs = XMLSimpleStr(npNLOval);
0264           sscanf(npNLOs, "%d", &npNLO);
0265         }
0266
0267         xmlEventNodes.resize(1);
0268         xmlEventNodes[0] = xmlEvent->getDocumentElement();
0269       }
0270       mode = kEvent;
0271     }
0272
0273     if (mode == kNone)
0274       throw cms::Exception("InvalidFormat") << "LHE file has invalid format" << std::endl;
0275
0276     buffer.clear();
0277   }
0278
0279   void LHEReader::XMLHandler::endElement(const XMLCh *const uri,
0280                                          const XMLCh *const localname,
0281                                          const XMLCh *const qname) {
0282     std::string name((const char *)XMLSimpleStr(qname));
0283
0284     if (mode) {
0285       if (mode == kHeader && xmlNodes.size() > 1) {
0286         xmlNodes.resize(xmlNodes.size() - 1);
0287         return;
0288       } else if (mode == kHeader) {
0289         std::unique_ptr<DOMLSSerializer> writer(impl->createLSSerializer());
0290         std::unique_ptr<DOMLSOutput> outputDesc(impl->createLSOutput());
0291         assert(outputDesc.get());
0292         outputDesc->setEncoding(XMLUniStr("UTF-8"));
0293
0294         for (DOMNode *node = xmlNodes[0]->getFirstChild(); node; node = node->getNextSibling()) {
0295           XMLSimpleStr buffer(writer->writeToString(node));
0296
0297           std::string type;
0298           const char *p, *q;
0299           DOMElement *elem;
0300
0301           switch (node->getNodeType()) {
0302             case DOMNode::ELEMENT_NODE:
0303               elem = static_cast<DOMElement *>(node);
0304               type = (const char *)XMLSimpleStr(elem->getTagName());
0305               p = std::strchr((const char *)buffer, '>') + 1;
0306               q = std::strrchr(p, '<');
0307               break;
0308             case DOMNode::COMMENT_NODE:
0309               type = "";
0310               p = buffer + 4;
0311               q = buffer + strlen(buffer) - 3;
0312               break;
0313             default:
0314               type = "<>";
0315               p = buffer + std::strspn(buffer, " \t\r\n");
0316               if (!*p)
0317                 continue;
0318               q = p + strlen(p);
0319           }
0320           LHERunInfo::Header header(type);
0321           fillHeader(header, p, q - p);
0322           headers.push_back(header);
0323         }
0324
0325         xmlHeader->release();
0326         xmlHeader = nullptr;
0327       } else if (name == "event" && mode == kEvent &&
0328                  (skipEvent || (!xmlEventNodes.empty()))) {  // handling of weights in LHE file
0329
0330         if (skipEvent) {
0331           gotObject = mode;
0332           mode = kNone;
0333           return;
0334         }
0335
0336         for (DOMNode *node = xmlEventNodes[0]->getFirstChild(); node; node = node->getNextSibling()) {
0337           switch (node->getNodeType()) {
0338             case DOMNode::ELEMENT_NODE:  // rwgt
0339               for (DOMNode *rwgt = xmlEventNodes[1]->getFirstChild(); rwgt; rwgt = rwgt->getNextSibling()) {
0340                 DOMNode *attr = rwgt->getAttributes()->item(0);
0341                 XMLSimpleStr atname(attr->getNodeValue());
0342                 XMLSimpleStr weight(rwgt->getFirstChild()->getNodeValue());
0343                 switch (rwgt->getNodeType()) {
0344                   case DOMNode::ELEMENT_NODE:
0345                     weightsinevent.push_back(std::make_pair((const char *)atname, (const char *)weight));
0346                     break;
0347                   default:
0348                     break;
0349                 }
0350               }
0351               break;
0352             case DOMNode::TEXT_NODE:  // event information
0353             {
0354               XMLSimpleStr data(node->getNodeValue());
0355               buffer.append(data);
0356             } break;
0357             default:
0358               break;
0359           }
0360         }
0361       } else if (mode == kEvent) {
0362         //skip unknown tags
0363         return;
0364       }
0365
0366       if (gotObject != kNone)
0367         throw cms::Exception("InvalidState") << "Unexpected pileup in"
0368                                                 " LHEReader::XMLHandler::endElement"
0369                                              << std::endl;
0370
0371       gotObject = mode;
0372       mode = kNone;
0373     }
0374   }
0375
0376   void LHEReader::XMLHandler::characters(const XMLCh *const data_, const XMLSize_t length) {
0377     if (mode == kHeader) {
0378       DOMText *text = xmlHeader->createTextNode(data_);
0379       xmlNodes.back()->appendChild(text);
0380       return;
0381     }
0382
0383     if (XMLSimpleStr::isAllSpaces(data_, length))
0384       return;
0385
0386     unsigned int offset = 0;
0387     while (offset < length && XMLSimpleStr::isSpace(data_[offset]))
0388       offset++;
0389
0390     if (mode == kEvent) {
0391       if (!skipEvent) {
0392         DOMText *text = xmlEvent->createTextNode(data_ + offset);
0393         xmlEventNodes.back()->appendChild(text);
0394       }
0395       return;
0396     }
0397
0398     if (mode == kNone)
0399       throw cms::Exception("InvalidFormat") << "LHE file has invalid format" << std::endl;
0400
0401     XMLSimpleStr data(data_ + offset);
0402     buffer.append(data);
0403   }
0404
0405   void LHEReader::XMLHandler::comment(const XMLCh *const data_, const XMLSize_t length) {
0406     if (mode == kHeader) {
0407       DOMComment *comment = xmlHeader->createComment(data_);
0408       xmlNodes.back()->appendChild(comment);
0409       return;
0410     }
0411
0412     XMLSimpleStr data(data_);
0413
0414     LHERunInfo::Header header;
0415     fillHeader(header, data);
0416     headers.push_back(header);
0417   }
0418
0419   LHEReader::LHEReader(const edm::ParameterSet &params)
0420       : fileURLs(params.getUntrackedParameter<std::vector<std::string> >("fileNames")),
0421         strName(""),
0422         firstEvent(params.getUntrackedParameter<unsigned int>("skipEvents", 0)),
0423         maxEvents(params.getUntrackedParameter<int>("limitEvents", -1)),
0424         curIndex(0),
0425         handler(new XMLHandler()) {}
0426
0427   LHEReader::LHEReader(const std::vector<std::string> &fileNames, unsigned int firstEvent)
0428       : fileURLs(fileNames),
0429         strName(""),
0430         firstEvent(firstEvent),
0431         maxEvents(-1),
0432         curIndex(0),
0433         handler(new XMLHandler()) {}
0434
0435   LHEReader::LHEReader(const std::string &inputs, unsigned int firstEvent)
0436       : strName(inputs), firstEvent(firstEvent), maxEvents(-1), curIndex(0), handler(new XMLHandler()) {}
0437
0438   LHEReader::~LHEReader() {
0439     // Explicitly release "orphaned" resources
0440     // that were created through DOM implementation
0441     // createXXXX factory method *before* last
0442     // XMLPlatformUtils::Terminate is called.
0443     handler.release();
0444     curDoc.release();
0445     curSource.release();
0446   }
0447
0448   std::shared_ptr<LHEEvent> LHEReader::next(bool *newFileOpened) {
0449     while (curDoc.get() || curIndex < fileURLs.size() || (fileURLs.empty() && !strName.empty())) {
0450       if (!curDoc.get()) {
0451         if (!platform) {
0452           //If we read multiple files, the XercesPlatform must live longer than any one
0453           // XMLDocument.
0454           platform = XMLDocument::platformHandle();
0455         }
0456         if (!fileURLs.empty()) {
0457           logFileAction("  Initiating request to open LHE file ", fileURLs[curIndex]);
0458           curSource = std::make_unique<FileSource>(fileURLs[curIndex]);
0459           logFileAction("  Successfully opened LHE file ", fileURLs[curIndex]);
0460           if (newFileOpened != nullptr)
0461             *newFileOpened = true;
0462           ++curIndex;
0463         } else if (!strName.empty()) {
0464           curSource = std::make_unique<StringSource>(strName);
0465         }
0466         handler->reset();
0467         curDoc.reset(curSource->createReader(*handler));
0468         curRunInfo.reset();
0469       }
0470       handler->skipEvent = firstEvent > 0;
0471
0472       XMLHandler::Object event = handler->gotObject;
0473       handler->gotObject = XMLHandler::kNone;
0474
0475       switch (event) {
0476         case XMLHandler::kNone:
0477           if (!curDoc->parse()) {
0478             curDoc.reset();
0479             logFileAction("  Closed LHE file ", fileURLs[curIndex - 1]);
0480             return std::shared_ptr<LHEEvent>();
0481           }
0482           break;
0483
0484         case XMLHandler::kHeader:
0485           break;
0486
0487         case XMLHandler::kInit: {
0488           std::istringstream data;
0489           data.str(handler->buffer);
0490           handler->buffer.clear();
0491
0492           curRunInfo = std::make_shared<LHERunInfo>(data);
0493
0494           std::for_each(handler->headers.begin(),
0495                         handler->headers.end(),
0496                         std::bind(&LHERunInfo::addHeader, curRunInfo.get(), std::placeholders::_1));
0497           handler->headers.clear();
0498         } break;
0499
0500         case XMLHandler::kComment:
0501           break;
0502
0503         case XMLHandler::kEvent: {
0504           if (!curRunInfo.get())
0505             throw cms::Exception("InvalidState") << "Got LHE event without"
0506                                                     " initialization."
0507                                                  << std::endl;
0508
0509           if (firstEvent > 0) {
0510             firstEvent--;
0511             continue;
0512           }
0513
0514           if (maxEvents == 0)
0515             return std::shared_ptr<LHEEvent>();
0516           else if (maxEvents > 0)
0517             maxEvents--;
0518
0519           std::istringstream data;
0520           data.str(handler->buffer);
0521           handler->buffer.clear();
0522
0523           std::shared_ptr<LHEEvent> lheevent;
0524           lheevent = std::make_shared<LHEEvent>(curRunInfo, data);
0525           const XMLHandler::wgt_info &info = handler->weightsinevent;
0526           for (size_t i = 0; i < info.size(); ++i) {
0527             double num = -1.0;
0528             sscanf(info[i].second.c_str(), "%le", &num);
0529             lheevent->addWeight(gen::WeightsInfo(info[i].first, num));
0530           }
0531           lheevent->setNpLO(handler->npLO);
0532           lheevent->setNpNLO(handler->npNLO);
0533           lheevent->setEvtNum(handler->evtnum);
0534           handler->evtnum = -1;
0535           //fill scales
0536           if (!handler->scales.empty()) {
0537             lheevent->setScales(handler->scales);
0538           }
0539           return lheevent;
0540         }
0541       }
0542     }
0543
0544     return std::shared_ptr<LHEEvent>();
0545   }
0546
0547 }  // namespace lhef