include/Pythia8Plugins/LHEH5v2.h

0001 // The LHEH5 code below has been adapted from
0002 // https://gitlab.com/hpcgen/tools authored by Holger Schulz
0003 // and Stefan Hoeche, and developed under the Scientific Discovery
0004 // through Advanced Computing (SciDAC) program funded by
0005 // the U.S. Department of Energy, Office of Science, Advanced Scientific
0006 // Computing Research.
0007 //
0008 // Note, this header can be used in conjuction with LHAHF5v2.h.
0009 //
0010 // Fermilab Software Legal Information (BSD License)
0011 // Copyright (c) 2009, FERMI NATIONAL ACCELERATOR LABORATORY
0012 // All rights reserved.
0013 //
0014 // Redistribution and use in source and binary forms, with or without
0015 // modification, are permitted provided that the following conditions
0016 // are met:
0017 //
0018 // Redistributions of source code must retain the above copyright
0019 // notice, this list of conditions and the following disclaimer.
0020 //
0021 // Redistributions in binary form must reproduce the above copyright
0022 // notice, this list of conditions and the following disclaimer in the
0023 // documentation and/or other materials provided with the
0024 // distribution.
0025 //
0026 // Neither the name of the FERMI NATIONAL ACCELERATOR LABORATORY, nor
0027 // the names of its contributors may be used to endorse or promote
0028 // products derived from this software without specific prior written
0029 // permission.
0030 //
0031 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
0032 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
0033 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
0034 // FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
0035 // COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
0036 // INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
0037 // (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
0038 // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
0039 // HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
0040 // STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
0041 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
0042 // OF THE POSSIBILITY OF SUCH DAMAGE.
0043
0044 #ifndef Pythia8_LHEH5v2_H
0045 #define Pythia8_LHEH5v2_H
0046
0047 // Standard includes.
0048 #include <iostream>
0049 #include <string>
0050 #include <vector>
0051
0052 // HighFive includes.
0053 #include "highfive/H5File.hpp"
0054 #include "highfive/H5FileDriver.hpp"
0055 #include "highfive/H5DataSet.hpp"
0056
0057 using namespace HighFive;
0058
0059 namespace LHEH5 {
0060
0061 //==========================================================================
0062
0063 // Process info struct.
0064
0065 //--------------------------------------------------------------------------
0066
0067 struct ProcInfo {
0068   int pid, nplo, npnlo;
0069   double unitwgt, xsec, error;
0070
0071   inline ProcInfo(int _pid,int _nplo,int _npnlo,
0072     double _unitwgt,double _xsec,double _error):
0073     pid(_pid), nplo(_nplo), npnlo(_npnlo),
0074       unitwgt(_unitwgt), xsec(_xsec), error(_error) {}
0075 };
0076
0077 //--------------------------------------------------------------------------
0078
0079 // Print process info.
0080
0081 std::ostream& operator<<(std::ostream& s,const ProcInfo& p) {
0082   return s<<"[pid="<<p.pid<<",nplo="<<p.nplo<<",npnlo="<<p.npnlo<<",xsec="
0083           <<p.xsec<<",err="<<p.error<<",unitwgt="<<p.unitwgt<<"]";
0084 }
0085
0086 //==========================================================================
0087
0088 // Particle struct.
0089
0090 //--------------------------------------------------------------------------
0091
0092 struct Particle {
0093   int id, st, mo1, mo2, cl1, cl2;
0094   double px, py, pz, e, m, lt, sp;
0095
0096   inline Particle(int _id,int _st,int _mo1,int _mo2,int _cl1,int _cl2,
0097     double _px,double _py,double _pz,double _e,double _m,
0098     double _lt,double _sp):
0099       id(_id), st(_st), mo1(_mo1), mo2(_mo2), cl1(_cl1), cl2(_cl2),
0100       px(_px), py(_py), pz(_pz), e(_e), m(_m), lt(_lt), sp(_sp) {}
0101   };// end of struct Particle
0102
0103 //--------------------------------------------------------------------------
0104
0105 // Print a particle
0106
0107 std::ostream& operator<<(std::ostream& s,const Particle& p) {
0108   return s << "{id=" << p.id << ",st=" << p.st
0109            << ",mo=[" << p.mo1 << "," << p.mo2 << "]"
0110            << ",cl=[" << p.cl1 << "," << p.cl2 << "]"
0111            << ",p=(" << p.e << "," << p.px << ","
0112            << p.py << "," << p.pz << ")}";
0113 }
0114
0115 //==========================================================================
0116
0117 // Event struct.
0118
0119 //--------------------------------------------------------------------------
0120
0121 struct Event : public std::vector<Particle> {
0122   ProcInfo pinfo;
0123   size_t trials;
0124   std::vector<double> wgts;
0125   double mur, muf, muq, aqed, aqcd;
0126   std::vector<Particle> ctparts;
0127   int ijt, kt, i, j, k;
0128   double z1, z2, bbpsw, psw;
0129
0130   inline Event(const ProcInfo &_pinfo,
0131     size_t _trials, std::vector<double> _wgts,
0132     double _mur, double _muf, double _muq,
0133     double _aqed, double _aqcd):
0134     pinfo(_pinfo), trials(_trials), wgts(_wgts),
0135       mur(_mur), muf(_muf), muq(_muq), aqed(_aqed), aqcd(_aqcd),
0136       ijt(-1), kt(-1), i(-1), j(-1), k(-1),
0137       z1(0), z2(0), bbpsw(0), psw(0) {}
0138   inline void AddCTInfo(int _ijt, int _kt, int _i, int _j, int _k,
0139     double _z1, double _z2, double _bbw, double _w) {
0140     ijt=_ijt; kt=_kt, i=_i; j=_j; k=_k;
0141     z1=_z1; z2=_z2; bbpsw=_bbw; psw=_w;
0142   }
0143
0144 };
0145
0146 //--------------------------------------------------------------------------
0147
0148 // Print an event.
0149
0150 std::ostream& operator<<(std::ostream& s, const Event& e) {
0151   s << "Event " << e.pinfo << " {\n"
0152     << "  trials=" << e.trials << ",weights=(" << e.wgts[0];
0153   for (size_t i(1); i<e.wgts.size(); ++i) s << "," << e.wgts[i];
0154   s << ")\n  mur=" << e.mur << ", muf=" << e.muf << ", muq=" << e.muq
0155     << ",aqed=" << e.aqed << ",aqcd=" << e.aqcd << "\n";
0156   for (size_t i(0); i<e.size(); ++i) s << "  " << e[i] << "\n";
0157   if (!e.ctparts.empty() || e.psw) {
0158     s << "  (" << e.ijt << "," << e.kt << ")->(" << e.i << "," << e.j
0159       << "," << e.k << "), z1=" << e.z1 << ", z2=" << e.z2
0160       << ", bbpsw=" << e.bbpsw << ", psw=" << e.psw << "\n";
0161     for (size_t i(0); i<e.ctparts.size(); ++i) s<<"  "<<e.ctparts[i]<<"\n";
0162   }
0163   return s<<"}";
0164 }
0165
0166 //==========================================================================
0167
0168 // LHEFile class.
0169
0170 //--------------------------------------------------------------------------
0171
0172 class LHEFile {
0173
0174  public:
0175
0176   inline const std::vector<int>& Version() const { return version; }
0177
0178   inline double TotalXS() const {
0179     double xs(0.);
0180     for (size_t i(0); i<pinfo.size(); ++i) xs += pinfo[i][3];
0181     return xs;
0182   }
0183
0184   inline double SumTrials() const {
0185     double trials(0.);
0186     for (size_t i(0); i<evts.size(); ++i) trials += evts[i][3];
0187     return trials;
0188   }
0189
0190   inline double UnitWeight() const {
0191     double wgt(0.);
0192     for (size_t i(0); i<pinfo.size(); ++i) wgt += pinfo[i][5];
0193     // For version 2.0.0 multiply by pb/GeV^2.
0194     if (version[0]==2 && version[1]==0 && version[2]==0) wgt*=3.89379656e8;
0195     return wgt;
0196   }
0197
0198   inline const std::vector<std::string>& WeightNames() const {
0199     return wgtnames;
0200   }
0201
0202   inline size_t NProcesses() const { return pinfo.size(); }
0203   inline ProcInfo GetProcInfo(const size_t pid) const {
0204     return ProcInfo(pid, pinfo[pid][1], pinfo[pid][2],
0205       pinfo[pid][5], pinfo[pid][3], pinfo[pid][4]);
0206   }
0207
0208   inline size_t NEvents() const { return evts.size(); }
0209   inline Event GetEvent(size_t i) const {
0210     std::vector<double> wgts(evts[i].begin()+9, evts[i].end());
0211     Event e(GetProcInfo(evts[i][0] ? evts[i][0]-1 : 0), evts[i][3], wgts,
0212       evts[i][6], evts[i][5], evts[i][4], evts[i][7], evts[i][8]);
0213     double wgt(0.);
0214     for (std::vector<double>::const_iterator it(wgts.begin());
0215          it!=wgts.end(); ++it) wgt += std::abs(*it);
0216     if (!wgt) return e;
0217     for (int n(0); n<evts[i][1]; ++n)
0218       e.push_back(GetParticle(evts[i][2]-evts[0][2]+n));
0219     if (!ctevts.empty()) {
0220       e.AddCTInfo(ctevts[i][0], ctevts[i][1], ctevts[i][2],
0221         ctevts[i][3], ctevts[i][4], ctevts[i][5],
0222         ctevts[i][6], ctevts[i][7], ctevts[i][8]);
0223       if (ctevts[i][0]>=0 && ctevts[i][1]>=0)
0224         for (int n(0); n<evts[i][1]+(ctevts[i][0]>=0 ? 1 : 0); ++n)
0225           e.ctparts.push_back(GetCTParticle(evts[i][2]-evts[0][2]+n));
0226     }
0227     return e;
0228   }
0229
0230   void Scale(const double s) {
0231     for (size_t i(0); i<evts.size(); ++i)
0232       for (size_t j(9); j<evts[i].size(); ++j) evts[i][j] *= s;
0233   }
0234
0235   void ReadHeader(File &file) {
0236     auto xfer_props = DataTransferProps{};
0237 #ifdef USING__MPI
0238     xfer_props.add(UseCollectiveIO{});
0239 #endif
0240     file.getDataSet("version").read(version,xfer_props);
0241     file.getDataSet("procInfo").read(pinfo,xfer_props);
0242     DataSet events(file.getDataSet("events"));
0243     auto attr_keys(events.listAttributeNames());
0244     Attribute a(events.getAttribute(attr_keys[0]));
0245     a.read(wgtnames);
0246     for (int i(0); i<9; ++i) wgtnames.erase(wgtnames.begin());
0247   }
0248
0249   void ReadEvents(File &file, size_t first_event, size_t n_events) {
0250     auto xfer_props = DataTransferProps{};
0251 #ifdef USING__MPI
0252     xfer_props.add(UseCollectiveIO{});
0253 #endif
0254     DataSet events(file.getDataSet("events"));
0255     std::vector<size_t> eoffsets{first_event, 0};
0256     std::vector<size_t> ecounts{n_events, 9+wgtnames.size()};
0257     evts.resize(n_events,std::vector<double>(9+wgtnames.size()));
0258     events.select(eoffsets,ecounts).read(evts, xfer_props);
0259     DataSet particles(file.getDataSet("particles"));
0260     std::vector<size_t> poffsets{(size_t)evts.front()[2], 0};
0261     size_t nmax(0);
0262     for (size_t i(0); i<pinfo.size(); ++i)
0263       nmax = std::max((size_t)std::max(pinfo[i][1],pinfo[i][2]+1), nmax);
0264     std::vector<size_t> pcounts{n_events*nmax, 13};
0265     parts.resize(n_events*nmax, std::vector<double>(13));
0266     particles.select(poffsets, pcounts).read(parts, xfer_props);
0267     if (file.exist("ctevents")) {
0268       DataSet ctevents(file.getDataSet("ctevents"));
0269       std::vector<size_t> cteoffsets{first_event, 0};
0270       std::vector<size_t> ctecounts{n_events, 9};
0271       ctevts.resize(n_events, std::vector<double>(9));
0272       ctevents.select(cteoffsets, ctecounts).read(ctevts, xfer_props);
0273       DataSet ctparticles(file.getDataSet("ctparticles"));
0274       std::vector<size_t> ctpoffsets{(size_t)evts.front()[2],0};
0275       std::vector<size_t> ctpcounts{n_events*nmax, 4};
0276       ctparts.resize(n_events*nmax, std::vector<double>(4));
0277       ctparticles.select(ctpoffsets, ctpcounts).read(ctparts, xfer_props);
0278     }
0279   }
0280
0281  private:
0282
0283   std::vector<int> version;
0284   std::vector<std::vector<double> > evts, parts, pinfo;
0285   std::vector<std::vector<double> > ctevts, ctparts;
0286   std::vector<std::string> wgtnames;
0287
0288   inline Particle GetParticle(size_t i) const {
0289     return Particle(parts[i][0], parts[i][1], parts[i][2], parts[i][3],
0290       parts[i][4], parts[i][5], parts[i][6], parts[i][7],
0291       parts[i][8], parts[i][9], parts[i][10],
0292       parts[i][11], parts[i][12]);
0293   }
0294
0295   inline Particle GetCTParticle(size_t i) const {
0296     return Particle(-1, -1, -1, -1, -1, -1, ctparts[i][0], ctparts[i][1],
0297       ctparts[i][2], ctparts[i][3], -1, -1, -1);
0298   }
0299
0300 };
0301
0302 //==========================================================================
0303
0304 }// end of namespace LHEH5
0305
0306 #endif // Pythia8_LHEH5v2_H