LHEReader.cc

Go to the documentation of this file.

#include <algorithm>
#include <iomanip>
#include <iostream>
#include <sstream>
#include <fstream>
#include <cstring>
#include <string>
#include <vector>
#include <cstdio>

#include <boost/bind.hpp>

#include <xercesc/sax2/Attributes.hpp>
#include <xercesc/dom/DOM.hpp>

#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Utilities/interface/Exception.h"
#include "FWCore/Utilities/interface/TimeOfDay.h"

#include "GeneratorInterface/LHEInterface/interface/LHEReader.h"
#include "GeneratorInterface/LHEInterface/interface/LHERunInfo.h"
#include "GeneratorInterface/LHEInterface/interface/LHEEvent.h"

#include "Utilities/StorageFactory/interface/IOTypes.h"
#include "Utilities/StorageFactory/interface/Storage.h"
#include "Utilities/StorageFactory/interface/StorageFactory.h"

#include "XMLUtils.h"

#include "boost/lexical_cast.hpp"

XERCES_CPP_NAMESPACE_USE

namespace lhef {

  static void logFileAction(char const* msg, std::string const& fileName) {
    edm::LogAbsolute("fileAction") << std::setprecision(0) << edm::TimeOfDay() << msg << fileName;
    edm::FlushMessageLog();
  }


class LHEReader::Source {
    public:
    Source() {}
    virtual ~Source() {}
    virtual XMLDocument *createReader(XMLDocument::Handler &handler) = 0;
};

class LHEReader::FileSource : public LHEReader::Source {
    public:
    FileSource(const std::string &fileURL)
    {
        auto storage =
            StorageFactory::get()->open(fileURL,
                                        IOFlags::OpenRead);

        if (!storage)
            throw cms::Exception("FileOpenError")
                << "Could not open LHE file \""
                << fileURL << "\" for reading"
                << std::endl;

        fileStream.reset(new StorageWrap(std::move(storage)));
    }

    ~FileSource() override {}

    XMLDocument *createReader(XMLDocument::Handler &handler) override
    { return new XMLDocument(fileStream, handler); }

    private:
    std::unique_ptr<StorageWrap>    fileStream;
};

class LHEReader::StringSource : public LHEReader::Source {
    public:
    StringSource(const std::string &inputs)
    {
        if (inputs.empty())
            throw cms::Exception("StreamOpenError")
                << "Empty LHE file string name \""
                << std::endl;

        std::stringstream * tmpis = new std::stringstream(inputs);
        fileStream.reset(tmpis);
    }

    ~StringSource() override {}

    XMLDocument *createReader(XMLDocument::Handler &handler) override
    { return new XMLDocument(fileStream, handler); }

    private:
  std::unique_ptr<std::istream> fileStream;
};

class LHEReader::XMLHandler : public XMLDocument::Handler {
    public:
  typedef std::vector<std::pair<std::string,std::string> > wgt_info;
    XMLHandler() :
        impl(nullptr),
        gotObject(kNone), mode(kNone),
        xmlHeader(nullptr), xmlEvent(nullptr), headerOk(false), npLO(-99), npNLO(-99) {}
    ~XMLHandler() override
    { if (xmlHeader) xmlHeader->release(); 
      if (xmlEvent) xmlEvent->release();   }

    enum Object {
        kNone = 0,
        kHeader,
        kInit,
        kComment,
        kEvent
    };

        void reset() { 
          headerOk = false; 
          weightsinevent.clear();
          gotObject = kNone;
          mode = kNone;
        }

        const wgt_info& weightInfo() const {return weightsinevent;}

    protected:
    void startElement(const XMLCh *const uri,
                      const XMLCh *const localname,
                      const XMLCh *const qname,
                      const Attributes &attributes) override;

    void endElement(const XMLCh *const uri,
                    const XMLCh *const localname,
                    const XMLCh *const qname) override;

        void characters (const XMLCh *const chars, const XMLSize_t length) override;
        void comment (const XMLCh *const chars, const XMLSize_t length) override;   

    private:
    friend class LHEReader;

        bool                            skipEvent = false;
        std::unique_ptr<DOMImplementation>  impl;
    std::string         buffer;
    Object              gotObject;
    Object              mode;
    DOMDocument         *xmlHeader;
        DOMDocument         *xmlEvent;
        std::vector<DOMElement*>    xmlNodes,xmlEventNodes;
    bool                headerOk;
    std::vector<LHERunInfo::Header> headers;
        wgt_info weightsinevent;
        int                             npLO;
        int                             npNLO;
        std::vector<float>              scales;
};

static void attributesToDom(DOMElement *dom, const Attributes &attributes)
{
    for(unsigned int i = 0; i < attributes.getLength(); i++) {
        const XMLCh *name = attributes.getQName(i);
        const XMLCh *value = attributes.getValue(i);

        dom->setAttribute(name, value);
    }
}

static void fillHeader(LHERunInfo::Header &header, const char *data,
                       int len = -1)
{
    const char *end = len >= 0 ? (data + len) : nullptr;
    while(*data && (!end || data < end)) {
        std::size_t len = std::strcspn(data, "\r\n");
        if (end && data + len > end)
            len = end - data;
        if (data[len] == '\r' && data[len + 1] == '\n')
            len += 2;
        else if (data[len])
            len++;
        header.addLine(std::string(data, len));
        data += len;
    }
}

void LHEReader::XMLHandler::startElement(const XMLCh *const uri,
                                         const XMLCh *const localname,
                                         const XMLCh *const qname,
                                         const Attributes &attributes)
{
  std::string name((const char*)XMLSimpleStr(qname));

  if (!headerOk) {
    if (name != "LesHouchesEvents")
      throw cms::Exception("InvalidFormat")
    << "LHE file has invalid header" << std::endl;
    headerOk = true;
    return;
  }
  
  if (mode == kHeader) {      
    DOMElement *elem = xmlHeader->createElement(qname);
    attributesToDom(elem, attributes);
    xmlNodes.back()->appendChild(elem);
    xmlNodes.push_back(elem);
    return;
  } else if ( mode == kEvent ) {

    if (skipEvent) {return;}

    DOMElement *elem = xmlEvent->createElement(qname);    
    attributesToDom(elem, attributes);

    //TODO this is a hack (even more than the rest of this class)
    if( name == "rwgt" ) {
      xmlEventNodes[0]->appendChild(elem);
    } else if (name == "wgt") {
      xmlEventNodes[1]->appendChild(elem);
    }
    else if (name == "scales") {
      for (XMLSize_t iscale=0; iscale<attributes.getLength(); ++iscale) {
        int ipart = 0;
        const char *scalename = XMLSimpleStr(attributes.getQName(iscale));
        int nmatch = sscanf(scalename,"pt_clust_%d",&ipart);
        
        if (nmatch!=1) {
          edm::LogError("Generator|LHEInterface")
          << "invalid attribute in <scales> tag"
          << std::endl;
        }

        float scaleval;
        const char *scalevalstr = XMLSimpleStr(attributes.getValue(iscale));
        sscanf(scalevalstr,"%e",&scaleval);
        
        scales.push_back(scaleval);
      }
    }
    xmlEventNodes.push_back(elem);
    return;
  } else if (mode == kInit) {
    //skip unknown tags in init block as well
    return;
  } else if (mode != kNone) {
    throw cms::Exception("InvalidFormat")
      << "LHE file has invalid format" << std::endl;
  }
  
  if (name == "header") {
    if (!impl)
      impl.reset(DOMImplementationRegistry::getDOMImplementation(XMLUniStr("Core")));

    xmlHeader = impl->createDocument(nullptr, qname, nullptr);
    xmlNodes.resize(1);
    xmlNodes[0] = xmlHeader->getDocumentElement();
    mode = kHeader;
  } if (name == "init") {
      mode = kInit;
  } else if (name == "event") {
    if (!skipEvent)
    {
      if (!impl)
    impl.reset(DOMImplementationRegistry::getDOMImplementation(XMLUniStr("Core")));

      if(xmlEvent)  xmlEvent->release();
      xmlEvent = impl->createDocument(nullptr, qname, nullptr);
      weightsinevent.resize(0);
      scales.clear();
    
      npLO = -99;
      npNLO = -99;
      const XMLCh *npLOval = attributes.getValue(XMLString::transcode("npLO"));
      if (npLOval) {
        const char *npLOs = XMLSimpleStr(npLOval);      
        sscanf(npLOs,"%d",&npLO);
      }
      const XMLCh *npNLOval = attributes.getValue(XMLString::transcode("npNLO"));
      if (npNLOval) {
        const char *npNLOs = XMLSimpleStr(npNLOval);      
         sscanf(npNLOs,"%d",&npNLO);
     }    
    
      xmlEventNodes.resize(1);
      xmlEventNodes[0] = xmlEvent->getDocumentElement();
    }
    mode = kEvent;    
  }
  
  if (mode == kNone)
    throw cms::Exception("InvalidFormat")
      << "LHE file has invalid format" << std::endl;
  
  buffer.clear();
}

void LHEReader::XMLHandler::endElement(const XMLCh *const uri,
                                       const XMLCh *const localname,
                                       const XMLCh *const qname)
{
  std::string name((const char*)XMLSimpleStr(qname));
  
  if (mode) {

    if (mode == kHeader && xmlNodes.size() > 1) {
      xmlNodes.resize(xmlNodes.size() - 1);
      return;
    } else if (mode == kHeader) {
      std::unique_ptr<DOMLSSerializer> writer(impl->createLSSerializer());
      std::unique_ptr<DOMLSOutput> outputDesc(impl->createLSOutput());
      assert(outputDesc.get());
      outputDesc->setEncoding(XMLUniStr("UTF-8"));

      for(DOMNode *node = xmlNodes[0]->getFirstChild();
      node; node = node->getNextSibling()) {
    XMLSimpleStr buffer(writer->writeToString(node));

    std::string type;
    const char *p, *q;
    DOMElement *elem;
    
    switch(node->getNodeType()) {
    case DOMNode::ELEMENT_NODE:   
        elem = static_cast<DOMElement*>(node);
        type = (const char*)XMLSimpleStr(
                         elem->getTagName());
        p = std::strchr((const char*)buffer,
                '>') + 1;
        q = std::strrchr(p, '<');    
      break;
    case DOMNode::COMMENT_NODE:
      type = "";
      p = buffer + 4;
      q = buffer + strlen(buffer) - 3;
      break;
    default:
      type = "<>";
      p = buffer +
        std::strspn(buffer, " \t\r\n");
      if (!*p)
        continue;
      q = p + strlen(p);
    }   
    LHERunInfo::Header header(type);
    fillHeader(header, p, q - p);
    headers.push_back(header);  
      }
      
      xmlHeader->release();
      xmlHeader = nullptr;
    }
    else if (name == "event" && 
    mode == kEvent && 
    (skipEvent || (!xmlEventNodes.empty()))) { // handling of weights in LHE file

      if (skipEvent)
      {
        gotObject = mode;
        mode = kNone;
        return;
      }

      for(DOMNode *node = xmlEventNodes[0]->getFirstChild();
      node; node = node->getNextSibling()) {
    switch( node->getNodeType() ) {
    case DOMNode::ELEMENT_NODE: // rwgt
      for(DOMNode *rwgt = xmlEventNodes[1]->getFirstChild();
          rwgt; rwgt = rwgt->getNextSibling()) {
        DOMNode* attr = rwgt->getAttributes()->item(0);
        XMLSimpleStr atname(attr->getNodeValue());
        XMLSimpleStr weight(rwgt->getFirstChild()->getNodeValue());
        switch( rwgt->getNodeType() ) {
        case DOMNode::ELEMENT_NODE:
          weightsinevent.push_back(std::make_pair((const char*)atname,
                              (const char*)weight));
          break;        
        default:              
          break;
        }     
      }
      break;
    case DOMNode::TEXT_NODE: // event information
      {
        XMLSimpleStr data(node->getNodeValue());
        buffer.append(data);
      }
      break;
    default:
      break;
    }
      }
    }
    else if (mode == kEvent) {
      //skip unknown tags
      return;
    }

    if (gotObject != kNone)
      throw cms::Exception("InvalidState")
    << "Unexpected pileup in"
    " LHEReader::XMLHandler::endElement"
    << std::endl;
    
    gotObject = mode;
    mode = kNone;
  }
}

void LHEReader::XMLHandler::characters(const XMLCh *const data_,
                                       const XMLSize_t length)
{
    if (mode == kHeader) {
        DOMText *text = xmlHeader->createTextNode(data_);
        xmlNodes.back()->appendChild(text);
        return;
    }   

    if (XMLSimpleStr::isAllSpaces(data_, length))
        return;

    unsigned int offset = 0;
    while(offset < length && XMLSimpleStr::isSpace(data_[offset]))
        offset++;   

    if( mode == kEvent ) {
          if (!skipEvent)
          {
            DOMText *text = xmlEvent->createTextNode(data_+offset);
            xmlEventNodes.back()->appendChild(text);
          }
      return;
    }   

    if (mode == kNone)
        throw cms::Exception("InvalidFormat")
            << "LHE file has invalid format" << std::endl;

    XMLSimpleStr data(data_ + offset);
    buffer.append(data);
}

void LHEReader::XMLHandler::comment(const XMLCh *const data_,
                                    const XMLSize_t length)
{
    if (mode == kHeader) {
        DOMComment *comment = xmlHeader->createComment(data_);
        xmlNodes.back()->appendChild(comment);
        return;
    }

    XMLSimpleStr data(data_);

    LHERunInfo::Header header;
    fillHeader(header, data);
    headers.push_back(header);
}

LHEReader::LHEReader(const edm::ParameterSet &params) :
    fileURLs(params.getUntrackedParameter< std::vector<std::string> >("fileNames")),
    strName(""),
    firstEvent(params.getUntrackedParameter<unsigned int>("skipEvents", 0)),
    maxEvents(params.getUntrackedParameter<int>("limitEvents", -1)),
    curIndex(0), handler(new XMLHandler())
{
}

LHEReader::LHEReader(const std::vector<std::string> &fileNames,
                     unsigned int firstEvent) :
  fileURLs(fileNames), strName(""), firstEvent(firstEvent), maxEvents(-1),
    curIndex(0), handler(new XMLHandler())
{
}

  LHEReader::LHEReader(const std::string &inputs,
                     unsigned int firstEvent) :
    strName(inputs), firstEvent(firstEvent), maxEvents(-1),
  curIndex(0), handler(new XMLHandler())
{
}

LHEReader::~LHEReader()
{
  // Explicitly release "orphaned" resources
  // that were created through DOM implementation
  // createXXXX factory method *before* last
  // XMLPlatformUtils::Terminate is called.
  handler.release();
  curDoc.release();
  curSource.release();
}

  std::shared_ptr<LHEEvent> LHEReader::next(bool* newFileOpened)
  {
    while(curDoc.get() || curIndex < fileURLs.size() || (fileURLs.empty() && !strName.empty() ) ) {
      if (!curDoc.get()) {
        if(!platform) {
          //If we read multiple files, the XercesPlatform must live longer than any one
          // XMLDocument.
          platform = XMLDocument::platformHandle();
        }
        if ( !fileURLs.empty() ) {
          logFileAction("  Initiating request to open LHE file ", fileURLs[curIndex]);
          curSource.reset(new FileSource(fileURLs[curIndex]));
          logFileAction("  Successfully opened LHE file ", fileURLs[curIndex]);
          if ( newFileOpened != nullptr ) *newFileOpened = true;
          ++curIndex;
        } else if ( !strName.empty() ) {
          curSource.reset(new StringSource(strName));
        }
        handler->reset();
        curDoc.reset(curSource->createReader(*handler));
        curRunInfo.reset();
      }
      handler->skipEvent = firstEvent > 0;
    
      XMLHandler::Object event = handler->gotObject;
      handler->gotObject = XMLHandler::kNone;
    
        
      switch(event) {
      case XMLHandler::kNone:
        if (!curDoc->parse()) {
          curDoc.reset();
          logFileAction("  Closed LHE file ", fileURLs[curIndex - 1]);
          return std::shared_ptr<LHEEvent>();
        }
        break;
        
      case XMLHandler::kHeader:
        break;
        
      case XMLHandler::kInit:
    {
    std::istringstream data;
    data.str(handler->buffer);
    handler->buffer.clear();

        curRunInfo.reset(new LHERunInfo(data));
        
        std::for_each(handler->headers.begin(),
                      handler->headers.end(),
                      boost::bind(&LHERunInfo::addHeader,
                                  curRunInfo.get(), _1));
        handler->headers.clear();
    }
        break;
        
      case XMLHandler::kComment:
        break;
        
      case XMLHandler::kEvent:
    {
        if (!curRunInfo.get())
          throw cms::Exception("InvalidState")
            << "Got LHE event without"
            " initialization." << std::endl;

        if (firstEvent > 0) {
          firstEvent--;
          continue;
        }
        
        if (maxEvents == 0)
          return std::shared_ptr<LHEEvent>();
        else if (maxEvents > 0)
          maxEvents--;

    std::istringstream data;
    data.str(handler->buffer);
    handler->buffer.clear();

    std::shared_ptr<LHEEvent> lheevent;
    lheevent.reset(new LHEEvent(curRunInfo, data));
    const XMLHandler::wgt_info& info = handler->weightsinevent;
    for( size_t i=0; i< info.size(); ++i ) {
      double num = -1.0;
      sscanf(info[i].second.c_str(),"%le",&num);      
      lheevent->addWeight(gen::WeightsInfo(info[i].first,num));
    }
    lheevent->setNpLO(handler->npLO);
        lheevent->setNpNLO(handler->npNLO);
        //fill scales
        if (!handler->scales.empty()) {
          lheevent->setScales(handler->scales);
        }
        return lheevent;
    }
      }
    }
    
    return std::shared_ptr<LHEEvent>();
  }
  
} // namespace lhef