File indexing completed on 2025-01-09 23:33:34
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
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;
0093 if (_eventsTotal % _eventsInBlock != 0)
0094 throw cms::Exception("ReadError") << "block size does not match HDF5 file" << std::endl;
0095 _blocksRead = 0;
0096
0097
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
0114 if (_blocksRead == 0 && maxEventsIn >= 0 && (unsigned long int)maxEventsIn < _eventsTotal)
0115 _eventsTotal = (unsigned long int)maxEventsIn;
0116
0117 if (firstEventIn > 0 && _blocksRead > 0)
0118 _eventsRead = firstEventIn - 1;
0119
0120 _blocksRead = _eventsRead / _eventsInBlock;
0121 }
0122
0123 void H5Handler::readBlock() {
0124
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
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 ¶ms)
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
0193 curDoc = true;
0194
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;
0217 std::vector<double> xSection;
0218 std::vector<double> error;
0219 std::vector<double> unitWeight;
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
0234 const HEPRUP heprup(tmprup);
0235
0236 curRunInfo = std::make_shared<LHERunInfo>(heprup);
0237
0238 }
0239
0240
0241 int npLO = curSource->handler->npLO;
0242 int npNLO = curSource->handler->npNLO;
0243
0244
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
0255 unsigned int ip = 0;
0256
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
0278 const HEPEUP hepeup(tmp);
0279
0280 lheevent = std::make_shared<LHEEvent>(curRunInfo, hepeup);
0281
0282
0283
0284
0285
0286
0287
0288
0289
0290 lheevent->setNpLO(npLO);
0291 lheevent->setNpNLO(npNLO);
0292
0293
0294
0295
0296 return lheevent;
0297 }
0298 return std::shared_ptr<LHEEvent>();
0299 }
0300
0301 }