include/Pythia8Plugins/LHAHDF5.h

0001 // LHAHDF5.h is a part of the PYTHIA event generator.
0002 // Copyright (C) 2024 Torbjorn Sjostrand.
0003 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
0004 // Please respect the MCnet Guidelines, see GUIDELINES for details.
0005 // Author: Christian Preuss, January 2021.
0006
0007 #ifndef Pythia8_LHAHDF5_H
0008 #define Pythia8_LHAHDF5_H
0009
0010 // Interface includes.
0011 #include "Pythia8Plugins/LHEH5.h"
0012
0013 // Generator includes.
0014 #include "Pythia8/Pythia.h"
0015
0016 namespace Pythia8 {
0017
0018 //==========================================================================
0019
0020 // HDF5 file reader. Converts to Pythia-internal events by acting as
0021 // replacement Les Houches Event reader.
0022
0023 class LHAupH5 : public Pythia8::LHAup {
0024
0025  public:
0026
0027   LHAupH5(HighFive::File* h5fileIn, size_t firstEventIn, size_t readSizeIn,
0028     string versionIn = "") :
0029     nTrialsSav(0), nRead(0), isOldFormat(h5fileIn->exist("/index")),
0030     particleSav(h5fileIn->getGroup("particle")),
0031     eventSav(h5fileIn->getGroup("event")),
0032     initSav(h5fileIn->getGroup("init")),
0033     procInfoSav(h5fileIn->getGroup("procInfo")),
0034     readSizeSav(readSizeIn), firstEventSav(firstEventIn), fileSav(h5fileIn) {
0035
0036     // Check if version number exists in event file.
0037     valid = false;
0038     if (h5fileIn->exist("/init/version")) {
0039       DataSet version = initSav.getDataSet("version");
0040       vector<int> versionNo;
0041       version.read(versionNo);
0042       versionSav = to_string(versionNo[0]) + "." + to_string(versionNo[1])
0043         + "." + to_string(versionNo[2]);
0044       isOldFormat = (versionSav=="0.1.0");
0045     } else if (isOldFormat) {
0046       // Based on the existence of the index group, we can guess this is 0.1.0.
0047       versionSav = "0.1.0";
0048     } else {
0049       // Have to rely on input if we could not guess the version.
0050       versionSav = versionIn;
0051       if (versionSav=="") {
0052         cout << " LHAupH5 error - cannot determine version number.\n";
0053         return;
0054       }
0055     }
0056     cout << " LHAupH5 version " << versionSav << ".\n";
0057
0058     // Check if version supports multiweights. (Starting from 1.0.0.)
0059     hasMultiWts = !(versionSav=="0.1.0" || versionSav=="0.2.0");
0060     cout << " LHAupH5 file format "
0061          << (hasMultiWts?"supports":"does not support")
0062          << " multi-weights." << endl;
0063
0064     hid_t dspace;
0065     if( !isOldFormat ) {
0066       DataSet npLO  = procInfoSav.getDataSet("npLO");
0067       DataSet npNLO = procInfoSav.getDataSet("npNLO");
0068       npLO.read(npLOSav);
0069       npNLO.read(npNLOSav);
0070       dspace = H5Dget_space(fileSav->getDataSet("event/start").getId());
0071     } else {
0072       indexSav      = fileSav->getGroup("index");
0073       dspace = H5Dget_space(fileSav->getDataSet("index/start").getId());
0074     }
0075     // Check if size of read is compatible.
0076     if (readSizeIn > H5Sget_simple_extent_npoints(dspace) ) {
0077       cout << "H5 size request incompatible with file.\n";
0078       return;
0079     }
0080
0081     // Check for multiweights.
0082     if ( fileSav->exist("event/weight") && hasMultiWts ) {
0083       DataSet weights = eventSav.getDataSet("weight");
0084       auto attr_keys = weights.listAttributeNames();
0085       Attribute a = weights.getAttribute(attr_keys[0]);
0086       a.read(weightsNames);
0087     }
0088
0089     // This reads and holds the information of readSize events,
0090     // starting from firstEvent.
0091     if (!isOldFormat) {
0092       lheEvts2Sav = LHEH5::readEvents2(particleSav, eventSav,
0093         firstEventSav, readSizeSav, npLOSav, npNLOSav, hasMultiWts);
0094
0095     } else lheEvtsSav = LHEH5::readEvents(indexSav, particleSav, eventSav,
0096         firstEventSav, readSizeSav);
0097     valid = true;
0098
0099   }
0100
0101   // Read and set the info from init and procInfo.
0102   bool setInit() override;
0103   bool setEvent(int idProc=0) override;
0104   void forceStrategy(int strategyIn) {setStrategy(strategyIn);}
0105   size_t nTrials() {return nTrialsSav;}
0106
0107  private:
0108
0109   // HDF5 file.
0110   HighFive::File* fileSav{};
0111
0112   // HDF5 Groups.
0113   HighFive::Group indexSav, particleSav, eventSav, initSav, procInfoSav;
0114
0115   // Events from HDF5 file.
0116   LHEH5::Events  lheEvtsSav;
0117   LHEH5::Events2 lheEvts2Sav;
0118
0119   // Info for reader.
0120   size_t         readSizeSav, firstEventSav, nTrialsSav;
0121   int            npLOSav, npNLOSav;
0122   bool           valid, isOldFormat, hasMultiWts;
0123   string         versionSav;
0124
0125   // Additional parameters.
0126   int nRead;
0127
0128   // Multiweight vector. Reset each event.
0129   vector<double> weightsSav;
0130   vector<string> weightsNames;
0131
0132   // Particle production scales.
0133   LHAscales scalesNow;
0134
0135 };
0136
0137 //--------------------------------------------------------------------------
0138
0139 // HDF5 file reader. Converts to Pythia-internal events by acting as
0140 // replacement Les Houches Event reader.
0141
0142 bool LHAupH5::setInit() {
0143   int beamA, beamB;
0144   double energyA, energyB;
0145   int PDFgroupA, PDFgroupB;
0146   int PDFsetA, PDFsetB;
0147
0148   if (!valid) return false;
0149   initSav.getDataSet("beamA").read(beamA);
0150   initSav.getDataSet("energyA").read(energyA);
0151   initSav.getDataSet("PDFgroupA").read(PDFgroupA);
0152   initSav.getDataSet("PDFsetA").read(PDFsetA);
0153
0154   initSav.getDataSet("beamB").read(beamB);
0155   initSav.getDataSet("energyB").read(energyB);
0156   initSav.getDataSet("PDFgroupB").read(PDFgroupB);
0157   initSav.getDataSet("PDFsetB").read(PDFsetB);
0158
0159   setBeamA(beamA, energyA, PDFgroupA, PDFsetA);
0160   setBeamB(beamB, energyB, PDFgroupB, PDFsetB);
0161
0162   int weightingStrategy = 3;
0163   initSav.getDataSet("weightingStrategy").read(weightingStrategy);
0164   setStrategy(weightingStrategy);
0165
0166   int nProcesses = 1;
0167   initSav.getDataSet("numProcesses").read(nProcesses);
0168
0169   vector<int> procId;
0170   vector<double> xSection, error, unitWeight;
0171   procInfoSav.getDataSet("procId").read(procId);
0172   procInfoSav.getDataSet("xSection").read(xSection);
0173   procInfoSav.getDataSet("error").read(error);
0174   procInfoSav.getDataSet("unitWeight").read(unitWeight);
0175   infoPtr->sigmaLHEFSave.resize(0);
0176   for (int iProc(0); iProc<nProcesses; ++iProc) {
0177     addProcess(procId[iProc], xSection[iProc], error[iProc],
0178                unitWeight[iProc]);
0179     infoPtr->sigmaLHEFSave.push_back(xSection[iProc]);
0180   }
0181   return true;
0182
0183 }
0184
0185 //--------------------------------------------------------------------------
0186
0187 // Read an event.
0188
0189 bool LHAupH5::setEvent(int idProc) {
0190
0191   // Equivalent of end of file.
0192   if (!valid) return false;
0193   if (nRead >= readSizeSav) return false;
0194   LHEH5::EventHeader evtHeader = !isOldFormat ?
0195     lheEvts2Sav.mkEventHeader(nRead) : lheEvtsSav.mkEventHeader(nRead);
0196   weightsSav = evtHeader.weights;
0197   nTrialsSav += evtHeader.trials;
0198   // Skip zero-weight events, but add trials.
0199   while (weightsSav[0] == 0. && nRead < readSizeSav - 1) {
0200     ++nRead;
0201     evtHeader  = !isOldFormat ?
0202       lheEvts2Sav.mkEventHeader(nRead) : lheEvtsSav.mkEventHeader(nRead);
0203     weightsSav = evtHeader.weights;
0204     nTrialsSav   += evtHeader.trials;
0205   }
0206   // Communicate event weight to Info.
0207   infoPtr->weightContainerPtr->setWeightNominal( weightsSav[0] );
0208   infoPtr->weightContainerPtr->weightsLHEF.bookVectors(
0209     weightsSav, weightsNames );
0210
0211   xwgtupSave = weightsSav[0];
0212   idprupSave = evtHeader.pid;
0213   scalupSave = evtHeader.scale;
0214   aqedupSave = evtHeader.aqed;
0215   aqcdupSave = evtHeader.aqcd;
0216   setProcess(idprupSave, xwgtupSave, scalupSave, aqedupSave, aqcdupSave);
0217   double scalein = -1.;
0218   vector<LHEH5::Particle> particles;
0219   particles = !isOldFormat ? lheEvts2Sav.mkEvent(nRead)
0220     : lheEvtsSav.mkEvent(nRead);
0221
0222   // Set particles.
0223   int nPtcls = 0;
0224   for (int iPtcl(0); iPtcl<int(particles.size()); ++iPtcl) {
0225     LHEH5::Particle ptcl = particles.at(iPtcl);
0226     if (ptcl.id == 0) continue;
0227     nPtcls++;
0228     if (iPtcl < 2) addParticle(ptcl.id, ptcl.status, 0, 0,
0229         ptcl.color1, ptcl.color2, ptcl.px, ptcl.py, ptcl.pz, ptcl.e, ptcl.m,
0230         ptcl.lifetime, ptcl.spin, scalein);
0231     else addParticle(ptcl.id, ptcl.status, ptcl.mother1, ptcl.mother2,
0232         ptcl.color1, ptcl.color2, ptcl.px, ptcl.py, ptcl.pz, ptcl.e, ptcl.m,
0233         ptcl.lifetime, ptcl.spin, scalein);
0234   }
0235   nupSave = nPtcls;
0236
0237   // Scale setting
0238   scalesNow.clear();
0239   scalesNow.muf   = evtHeader.fscale;
0240   scalesNow.mur   = evtHeader.rscale;
0241   scalesNow.mups  = evtHeader.scale;
0242   infoPtr->scales = &scalesNow;
0243   ++nRead;
0244   return true;
0245
0246 }
0247
0248 //==========================================================================
0249
0250 } // end namespace Pythia8
0251
0252 #endif // Pythia8_LHAHDF5_H