LHEInterface/src/LH5Reader.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 "FWCore/MessageLogger/interface/MessageLogger.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/Utilities/interface/Exception.h"
0016 #include "FWCore/Utilities/interface/TimeOfDay.h"
0017
0018 #include "GeneratorInterface/LHEInterface/interface/LH5Reader.h"
0019 #include "GeneratorInterface/LHEInterface/interface/LHERunInfo.h"
0020 #include "GeneratorInterface/LHEInterface/interface/LHEEvent.h"
0021 //#include "GeneratorInterface/LHEInterface/code/HighFive/include/highfive/H5File.hpp"
0022 #include "GeneratorInterface/LHEInterface/interface/lheh5.h"
0023
0024 #include "Utilities/StorageFactory/interface/IOTypes.h"
0025 #include "Utilities/StorageFactory/interface/Storage.h"
0026 #include "Utilities/StorageFactory/interface/StorageFactory.h"
0027
0028 namespace lhef {
0029   using namespace lheh5;
0030
0031   static void logFileAction(char const *msg, std::string const &fileName) {
0032     edm::LogAbsolute("fileAction") << std::setprecision(0) << edm::TimeOfDay() << msg << fileName;
0033     edm::FlushMessageLog();
0034   }
0035
0036   class LH5Reader::Source {
0037   public:
0038     Source() {}
0039     virtual ~Source() {}
0040     std::unique_ptr<H5Handler> handler;
0041
0042   private:
0043   };
0044
0045   class LH5Reader::FileSource : public LH5Reader::Source {
0046   public:
0047     FileSource(const std::string &fileURL) {
0048       const std::string fileName = fileURL.substr(5, fileURL.length());
0049
0050       H5Handler *tmph = new H5Handler(fileName);
0051       handler.reset(tmph);
0052     }
0053
0054     ~FileSource() override {}
0055   };
0056
0057   class LH5Reader::StringSource : public LH5Reader::Source {
0058   public:
0059     StringSource(const std::string &inputs) {
0060       if (inputs.empty())
0061         throw cms::Exception("StreamOpenError") << "Empty LHE file string name \"" << std::endl;
0062
0063       H5Handler *tmph = new H5Handler(inputs);
0064       handler.reset(tmph);
0065     }
0066
0067     ~StringSource() override {}
0068   };
0069
0070   H5Handler::H5Handler(const std::string &fileNameIn)
0071       : h5file(new HighFive::File(fileNameIn)),
0072         indexStatus(h5file->exist("/index")),
0073         _index(h5file->getGroup(indexStatus ? "index" : "event")),
0074         _particle(h5file->getGroup("particle")),
0075         _event(h5file->getGroup("event")),
0076         _init(h5file->getGroup("init")),
0077         _procInfo(h5file->getGroup("procInfo")) {
0078     hid_t dspace;
0079     _formatType = 1;
0080
0081     if (indexStatus) {
0082       dspace = H5Dget_space(h5file->getDataSet("index/start").getId());
0083       _formatType = 2;
0084     } else {
0085       _index = h5file->getGroup("event");
0086       dspace = H5Dget_space(h5file->getDataSet("event/start").getId());
0087       _formatType = 1;
0088     }
0089
0090     _eventsTotal = H5Sget_simple_extent_npoints(dspace);
0091     _eventsRead = 0;
0092     _eventsInBlock = 6400 / 4;  // configurable parameter??
0093     if (_eventsTotal % _eventsInBlock != 0)
0094       throw cms::Exception("ReadError") << "block size does not match HDF5 file" << std::endl;
0095     _blocksRead = 0;
0096
0097     // Check if the file contains the npLO, npNLO information
0098     bool status = h5file->exist("/procInfo/npLO");
0099     if (status) {
0100       HighFive::DataSet _npLO = _procInfo.getDataSet("npLO");
0101       _npLO.read(npLO);
0102       HighFive::DataSet _npNLO = _procInfo.getDataSet("npNLO");
0103       _npNLO.read(npNLO);
0104     } else {
0105       npLO = 1;
0106       npNLO = 0;
0107     }
0108   };
0109
0110   void H5Handler::counter(const int firstEventIn, const int maxEventsIn) {
0111     if (maxEventsIn > 0 && firstEventIn > maxEventsIn)
0112       throw cms::Exception("ConfigurationError") << "\" firstEvent > maxEvents \"" << std::endl;
0113     // Reset maximum number of events to read
0114     if (_blocksRead == 0 && maxEventsIn >= 0 && (unsigned long int)maxEventsIn < _eventsTotal)
0115       _eventsTotal = (unsigned long int)maxEventsIn;
0116     // If there are multiple files, assume you only want to jump through the first file
0117     if (firstEventIn > 0 && _blocksRead > 0)
0118       _eventsRead = firstEventIn - 1;
0119     // Set blocks read to be in the correct place
0120     _blocksRead = _eventsRead / _eventsInBlock;
0121   }
0122
0123   void H5Handler::readBlock() {
0124     // Start counting at 0
0125     size_t iStart = _blocksRead * _eventsInBlock;
0126     size_t nEvents = _eventsInBlock;
0127     if ((unsigned long int)(iStart + nEvents) > _eventsTotal)
0128       return;
0129     _events1 = readEvents(_index, _particle, _event, iStart, nEvents);
0130     _blocksRead++;
0131   }
0132
0133   std::vector<lheh5::Particle> H5Handler::getEvent() {
0134     std::vector<lheh5::Particle> _vE;
0135     if (_eventsRead > _eventsTotal - 1)
0136       return std::vector<lheh5::Particle>();
0137     int checkEvents = (_blocksRead)*_eventsInBlock - 1;
0138     _vE = _events1.mkEvent(_eventsRead);
0139     _eventsRead++;
0140     if (checkEvents >= 0 && _eventsRead > (unsigned long int)checkEvents)
0141       readBlock();
0142     return _vE;
0143   }
0144
0145   lheh5::EventHeader H5Handler::getHeader() {
0146     // fragile, must be called before getEvent
0147     return _events1.mkEventHeader(_eventsRead);
0148   }
0149
0150   std::pair<lheh5::EventHeader, std::vector<lheh5::Particle> > H5Handler::getEventProperties() {
0151     std::vector<lheh5::Particle> _vE;
0152     lheh5::EventHeader _h;
0153     _h.nparticles = -1;
0154     if (_eventsRead > _eventsTotal - 1)
0155       return std::make_pair(_h, _vE);
0156     int checkEvents = (_blocksRead)*_eventsInBlock - 1;
0157     _vE = _events1.mkEvent(_eventsRead % _eventsInBlock);
0158     _h = _events1.mkEventHeader(_eventsRead % _eventsInBlock);
0159     _eventsRead++;
0160     if (checkEvents >= 0 && _eventsRead > (unsigned long int)checkEvents)
0161       readBlock();
0162     return std::make_pair(_h, _vE);
0163   }
0164
0165   LH5Reader::LH5Reader(const edm::ParameterSet &params)
0166       : fileURLs(params.getUntrackedParameter<std::vector<std::string> >("fileNames")),
0167         strName(""),
0168         firstEvent(params.getUntrackedParameter<unsigned int>("skipEvents", 0)),
0169         maxEvents(params.getUntrackedParameter<int>("limitEvents", -1)) {}
0170
0171   LH5Reader::LH5Reader(const std::vector<std::string> &fileNames, unsigned int firstEvent, int maxEvents)
0172       : fileURLs(fileNames), strName(""), firstEvent(firstEvent), maxEvents(maxEvents), curIndex(0), curDoc(false) {}
0173
0174   LH5Reader::LH5Reader(const std::string &inputs, unsigned int firstEvent, int maxEvents)
0175       : strName(inputs), firstEvent(firstEvent), maxEvents(maxEvents), curIndex(0) {}
0176
0177   LH5Reader::~LH5Reader() { curSource.release(); }
0178
0179   std::shared_ptr<LHEEvent> LH5Reader::next(bool *newFileOpened) {
0180     while (curDoc || curIndex < fileURLs.size() || (fileURLs.empty() && !strName.empty())) {
0181       if (!curDoc) {
0182         if (!fileURLs.empty()) {
0183           logFileAction("  Initiating request to open LHE file ", fileURLs[curIndex]);
0184           curSource = std::make_unique<FileSource>(fileURLs[curIndex]);
0185           logFileAction("  Successfully opened LHE file ", fileURLs[curIndex]);
0186           if (newFileOpened != nullptr)
0187             *newFileOpened = true;
0188           ++curIndex;
0189         } else if (!strName.empty()) {
0190           curSource = std::make_unique<StringSource>(strName);
0191         }
0192         // New "doc" has been opened.    This is the same as a new source.
0193         curDoc = true;
0194         // Set maxEvents and firstEvent
0195         curSource->handler->counter(firstEvent, maxEvents);
0196         curSource->handler->readBlock();
0197
0198         curRunInfo.reset();
0199         HEPRUP tmprup;
0200         int beamA, beamB;
0201         curSource->handler->_init.getDataSet("beamA").read(beamA);
0202         curSource->handler->_init.getDataSet("beamB").read(beamB);
0203         tmprup.IDBMUP = std::make_pair(beamA, beamB);
0204         double energyA, energyB;
0205         curSource->handler->_init.getDataSet("energyA").read(energyA);
0206         curSource->handler->_init.getDataSet("energyB").read(energyB);
0207         tmprup.EBMUP = std::make_pair(energyA, energyB);
0208         int PDFsetA, PDFsetB;
0209         curSource->handler->_init.getDataSet("PDFsetA").read(PDFsetA);
0210         curSource->handler->_init.getDataSet("PDFsetB").read(PDFsetB);
0211         tmprup.PDFSUP = std::make_pair(PDFsetA, PDFsetB);
0212         int PDFgroupA, PDFgroupB;
0213         curSource->handler->_init.getDataSet("PDFgroupA").read(PDFgroupA);
0214         curSource->handler->_init.getDataSet("PDFgroupB").read(PDFgroupB);
0215         tmprup.PDFGUP = std::make_pair(PDFgroupA, PDFgroupB);
0216         std::vector<int> procId;         // NOTE: C++17 allows int[numProcesses]
0217         std::vector<double> xSection;    // NOTE: C++17 allows double[numProcesses]
0218         std::vector<double> error;       // NOTE: C++17 allows double[numProcesses]
0219         std::vector<double> unitWeight;  // NOTE: C++17 allows double[numProcesses]
0220
0221         curSource->handler->_procInfo.getDataSet("procId").read(procId);
0222         curSource->handler->_procInfo.getDataSet("xSection").read(xSection);
0223         curSource->handler->_procInfo.getDataSet("error").read(error);
0224         curSource->handler->_procInfo.getDataSet("unitWeight").read(unitWeight);
0225
0226         tmprup.LPRUP = procId;
0227         tmprup.XSECUP = xSection;
0228         tmprup.XERRUP = error;
0229         tmprup.XMAXUP = unitWeight;
0230         tmprup.IDWTUP = 3;
0231         size_t numProcesses = procId.size();
0232         tmprup.NPRUP = numProcesses;
0233         // Use temporary process info block to define const HEPRUP
0234         const HEPRUP heprup(tmprup);
0235
0236         curRunInfo.reset(new LHERunInfo(heprup));
0237         // Run info has now been set when a new file is encountered
0238       }
0239       // Handler should be modified to have these capabilities
0240       // Maybe this is set event by event??
0241       int npLO = curSource->handler->npLO;
0242       int npNLO = curSource->handler->npNLO;
0243
0244       // Event-loop here
0245       std::pair<EventHeader, std::vector<Particle> > evp = curSource->handler->getEventProperties();
0246       EventHeader hd = evp.first;
0247       if (hd.nparticles < 0) {
0248         curDoc = false;
0249         logFileAction("  Closed LHE file ", fileURLs[curIndex - 1]);
0250         return std::shared_ptr<LHEEvent>();
0251       }
0252       HEPEUP tmp;
0253       tmp.resize(hd.nparticles);
0254       //Particle loop
0255       unsigned int ip = 0;
0256       //      for (auto part: curSource->handler->_events1.mkEvent(i)) {
0257       for (auto part : evp.second) {
0258         tmp.IDUP[ip] = part.id;
0259         tmp.ISTUP[ip] = part.status;
0260         tmp.MOTHUP[ip] = std::make_pair(part.mother1, part.mother2);
0261         tmp.ICOLUP[ip] = std::make_pair(part.color1, part.color2);
0262         tmp.VTIMUP[ip] = part.lifetime;
0263         tmp.SPINUP[ip] = part.spin;
0264         tmp.PUP[ip][0] = part.px;
0265         tmp.PUP[ip][1] = part.py;
0266         tmp.PUP[ip][2] = part.pz;
0267         tmp.PUP[ip][3] = part.e;
0268         tmp.PUP[ip][4] = part.m;
0269         ip++;
0270       }
0271       tmp.IDPRUP = hd.pid;
0272       tmp.XWGTUP = hd.weight;
0273       tmp.SCALUP = hd.scale;
0274       tmp.AQEDUP = hd.aqed;
0275       tmp.AQCDUP = hd.aqcd;
0276       std::shared_ptr<LHEEvent> lheevent;
0277       // Use temporary event to construct const HEPEUP;
0278       const HEPEUP hepeup(tmp);
0279
0280       lheevent.reset(new LHEEvent(curRunInfo, hepeup));
0281       // Might have to add this capability later
0282       /*          const XMLHandler::wgt_info &info = handler->weightsinevent;
0283           for (size_t i = 0; i < info.size(); ++i) {
0284             double num = -1.0;
0285             sscanf(info[i].second.c_str(), "%le", &num);
0286             lheevent->addWeight(gen::WeightsInfo(info[i].first, num));
0287             }*/
0288       // Currently these are set just at the beginning?
0289       // might be an event property
0290       lheevent->setNpLO(npLO);
0291       lheevent->setNpNLO(npNLO);
0292       //fill scales
0293       /*          if (!handler->scales.empty()) {
0294             lheevent->setScales(handler->scales);
0295             }*/
0296       return lheevent;
0297     }
0298     return std::shared_ptr<LHEEvent>();
0299   }
0300
0301 }  // namespace lhef