PYTHIA  8.312
LHEH5v2.h
1 // The LHEH5 code below has been adapted from
2 // https://gitlab.com/hpcgen/tools authored by Holger Schulz
3 // and Stefan Hoeche, and developed under the Scientific Discovery
4 // through Advanced Computing (SciDAC) program funded by
5 // the U.S. Department of Energy, Office of Science, Advanced Scientific
6 // Computing Research.
7 //
8 // Note, this header can be used in conjuction with LHAHF5v2.h.
9 //
10 // Fermilab Software Legal Information (BSD License)
11 // Copyright (c) 2009, FERMI NATIONAL ACCELERATOR LABORATORY
12 // All rights reserved.
13 //
14 // Redistribution and use in source and binary forms, with or without
15 // modification, are permitted provided that the following conditions
16 // are met:
17 //
18 // Redistributions of source code must retain the above copyright
19 // notice, this list of conditions and the following disclaimer.
20 //
21 // Redistributions in binary form must reproduce the above copyright
22 // notice, this list of conditions and the following disclaimer in the
23 // documentation and/or other materials provided with the
24 // distribution.
25 //
26 // Neither the name of the FERMI NATIONAL ACCELERATOR LABORATORY, nor
27 // the names of its contributors may be used to endorse or promote
28 // products derived from this software without specific prior written
29 // permission.
30 //
31 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
32 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
33 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
34 // FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
35 // COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
36 // INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
37 // (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
38 // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
39 // HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
40 // STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
41 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
42 // OF THE POSSIBILITY OF SUCH DAMAGE.
43 
44 #ifndef Pythia8_LHEH5v2_H
45 #define Pythia8_LHEH5v2_H
46 
47 // Standard includes.
48 #include <iostream>
49 #include <string>
50 #include <vector>
51 
52 // HighFive includes.
53 #include "highfive/H5File.hpp"
54 #include "highfive/H5FileDriver.hpp"
55 #include "highfive/H5DataSet.hpp"
56 
57 using namespace HighFive;
58 
59 namespace LHEH5 {
60 
61 //==========================================================================
62 
63 // Process info struct.
64 
65 //--------------------------------------------------------------------------
66 
67 struct ProcInfo {
68  int pid, nplo, npnlo;
69  double unitwgt, xsec, error;
70 
71  inline ProcInfo(int _pid,int _nplo,int _npnlo,
72  double _unitwgt,double _xsec,double _error):
73  pid(_pid), nplo(_nplo), npnlo(_npnlo),
74  unitwgt(_unitwgt), xsec(_xsec), error(_error) {}
75 };
76 
77 //--------------------------------------------------------------------------
78 
79 // Print process info.
80 
81 std::ostream& operator<<(std::ostream& s,const ProcInfo& p) {
82  return s<<"[pid="<<p.pid<<",nplo="<<p.nplo<<",npnlo="<<p.npnlo<<",xsec="
83  <<p.xsec<<",err="<<p.error<<",unitwgt="<<p.unitwgt<<"]";
84 }
85 
86 //==========================================================================
87 
88 // Particle struct.
89 
90 //--------------------------------------------------------------------------
91 
92 struct Particle {
93  int id, st, mo1, mo2, cl1, cl2;
94  double px, py, pz, e, m, lt, sp;
95 
96  inline Particle(int _id,int _st,int _mo1,int _mo2,int _cl1,int _cl2,
97  double _px,double _py,double _pz,double _e,double _m,
98  double _lt,double _sp):
99  id(_id), st(_st), mo1(_mo1), mo2(_mo2), cl1(_cl1), cl2(_cl2),
100  px(_px), py(_py), pz(_pz), e(_e), m(_m), lt(_lt), sp(_sp) {}
101  };// end of struct Particle
102 
103 //--------------------------------------------------------------------------
104 
105 // Print a particle
106 
107 std::ostream& operator<<(std::ostream& s,const Particle& p) {
108  return s << "{id=" << p.id << ",st=" << p.st
109  << ",mo=[" << p.mo1 << "," << p.mo2 << "]"
110  << ",cl=[" << p.cl1 << "," << p.cl2 << "]"
111  << ",p=(" << p.e << "," << p.px << ","
112  << p.py << "," << p.pz << ")}";
113 }
114 
115 //==========================================================================
116 
117 // Event struct.
118 
119 //--------------------------------------------------------------------------
120 
121 struct Event : public std::vector<Particle> {
122  ProcInfo pinfo;
123  size_t trials;
124  std::vector<double> wgts;
125  double mur, muf, muq, aqed, aqcd;
126  std::vector<Particle> ctparts;
127  int ijt, kt, i, j, k;
128  double z1, z2, bbpsw, psw;
129 
130  inline Event(const ProcInfo &_pinfo,
131  size_t _trials, std::vector<double> _wgts,
132  double _mur, double _muf, double _muq,
133  double _aqed, double _aqcd):
134  pinfo(_pinfo), trials(_trials), wgts(_wgts),
135  mur(_mur), muf(_muf), muq(_muq), aqed(_aqed), aqcd(_aqcd),
136  ijt(-1), kt(-1), i(-1), j(-1), k(-1),
137  z1(0), z2(0), bbpsw(0), psw(0) {}
138  inline void AddCTInfo(int _ijt, int _kt, int _i, int _j, int _k,
139  double _z1, double _z2, double _bbw, double _w) {
140  ijt=_ijt; kt=_kt, i=_i; j=_j; k=_k;
141  z1=_z1; z2=_z2; bbpsw=_bbw; psw=_w;
142  }
143 
144 };
145 
146 //--------------------------------------------------------------------------
147 
148 // Print an event.
149 
150 std::ostream& operator<<(std::ostream& s, const Event& e) {
151  s << "Event " << e.pinfo << " {\n"
152  << " trials=" << e.trials << ",weights=(" << e.wgts[0];
153  for (size_t i(1); i<e.wgts.size(); ++i) s << "," << e.wgts[i];
154  s << ")\n mur=" << e.mur << ", muf=" << e.muf << ", muq=" << e.muq
155  << ",aqed=" << e.aqed << ",aqcd=" << e.aqcd << "\n";
156  for (size_t i(0); i<e.size(); ++i) s << " " << e[i] << "\n";
157  if (!e.ctparts.empty() || e.psw) {
158  s << " (" << e.ijt << "," << e.kt << ")->(" << e.i << "," << e.j
159  << "," << e.k << "), z1=" << e.z1 << ", z2=" << e.z2
160  << ", bbpsw=" << e.bbpsw << ", psw=" << e.psw << "\n";
161  for (size_t i(0); i<e.ctparts.size(); ++i) s<<" "<<e.ctparts[i]<<"\n";
162  }
163  return s<<"}";
164 }
165 
166 //==========================================================================
167 
168 // LHEFile class.
169 
170 //--------------------------------------------------------------------------
171 
172 class LHEFile {
173 
174  public:
175 
176  inline const std::vector<int>& Version() const { return version; }
177 
178  inline double TotalXS() const {
179  double xs(0.);
180  for (size_t i(0); i<pinfo.size(); ++i) xs += pinfo[i][3];
181  return xs;
182  }
183 
184  inline double SumTrials() const {
185  double trials(0.);
186  for (size_t i(0); i<evts.size(); ++i) trials += evts[i][3];
187  return trials;
188  }
189 
190  inline double UnitWeight() const {
191  double wgt(0.);
192  for (size_t i(0); i<pinfo.size(); ++i) wgt += pinfo[i][5];
193  // For version 2.0.0 multiply by pb/GeV^2.
194  if (version[0]==2 && version[1]==0 && version[2]==0) wgt*=3.89379656e8;
195  return wgt;
196  }
197 
198  inline const std::vector<std::string>& WeightNames() const {
199  return wgtnames;
200  }
201 
202  inline size_t NProcesses() const { return pinfo.size(); }
203  inline ProcInfo GetProcInfo(const size_t pid) const {
204  return ProcInfo(pid, pinfo[pid][1], pinfo[pid][2],
205  pinfo[pid][5], pinfo[pid][3], pinfo[pid][4]);
206  }
207 
208  inline size_t NEvents() const { return evts.size(); }
209  inline Event GetEvent(size_t i) const {
210  std::vector<double> wgts(evts[i].begin()+9, evts[i].end());
211  Event e(GetProcInfo(evts[i][0] ? evts[i][0]-1 : 0), evts[i][3], wgts,
212  evts[i][6], evts[i][5], evts[i][4], evts[i][7], evts[i][8]);
213  double wgt(0.);
214  for (std::vector<double>::const_iterator it(wgts.begin());
215  it!=wgts.end(); ++it) wgt += std::abs(*it);
216  if (!wgt) return e;
217  for (int n(0); n<evts[i][1]; ++n)
218  e.push_back(GetParticle(evts[i][2]-evts[0][2]+n));
219  if (!ctevts.empty()) {
220  e.AddCTInfo(ctevts[i][0], ctevts[i][1], ctevts[i][2],
221  ctevts[i][3], ctevts[i][4], ctevts[i][5],
222  ctevts[i][6], ctevts[i][7], ctevts[i][8]);
223  if (ctevts[i][0]>=0 && ctevts[i][1]>=0)
224  for (int n(0); n<evts[i][1]+(ctevts[i][0]>=0 ? 1 : 0); ++n)
225  e.ctparts.push_back(GetCTParticle(evts[i][2]-evts[0][2]+n));
226  }
227  return e;
228  }
229 
230  void Scale(const double s) {
231  for (size_t i(0); i<evts.size(); ++i)
232  for (size_t j(9); j<evts[i].size(); ++j) evts[i][j] *= s;
233  }
234 
235  void ReadHeader(File &file) {
236  auto xfer_props = DataTransferProps{};
237 #ifdef USING__MPI
238  xfer_props.add(UseCollectiveIO{});
239 #endif
240  file.getDataSet("version").read(version,xfer_props);
241  file.getDataSet("procInfo").read(pinfo,xfer_props);
242  DataSet events(file.getDataSet("events"));
243  auto attr_keys(events.listAttributeNames());
244  Attribute a(events.getAttribute(attr_keys[0]));
245  a.read(wgtnames);
246  for (int i(0); i<9; ++i) wgtnames.erase(wgtnames.begin());
247  }
248 
249  void ReadEvents(File &file, size_t first_event, size_t n_events) {
250  auto xfer_props = DataTransferProps{};
251 #ifdef USING__MPI
252  xfer_props.add(UseCollectiveIO{});
253 #endif
254  DataSet events(file.getDataSet("events"));
255  std::vector<size_t> eoffsets{first_event, 0};
256  std::vector<size_t> ecounts{n_events, 9+wgtnames.size()};
257  evts.resize(n_events,std::vector<double>(9+wgtnames.size()));
258  events.select(eoffsets,ecounts).read(evts, xfer_props);
259  DataSet particles(file.getDataSet("particles"));
260  std::vector<size_t> poffsets{(size_t)evts.front()[2], 0};
261  size_t nmax(0);
262  for (size_t i(0); i<pinfo.size(); ++i)
263  nmax = std::max((size_t)std::max(pinfo[i][1],pinfo[i][2]+1), nmax);
264  std::vector<size_t> pcounts{n_events*nmax, 13};
265  parts.resize(n_events*nmax, std::vector<double>(13));
266  particles.select(poffsets, pcounts).read(parts, xfer_props);
267  if (file.exist("ctevents")) {
268  DataSet ctevents(file.getDataSet("ctevents"));
269  std::vector<size_t> cteoffsets{first_event, 0};
270  std::vector<size_t> ctecounts{n_events, 9};
271  ctevts.resize(n_events, std::vector<double>(9));
272  ctevents.select(cteoffsets, ctecounts).read(ctevts, xfer_props);
273  DataSet ctparticles(file.getDataSet("ctparticles"));
274  std::vector<size_t> ctpoffsets{(size_t)evts.front()[2],0};
275  std::vector<size_t> ctpcounts{n_events*nmax, 4};
276  ctparts.resize(n_events*nmax, std::vector<double>(4));
277  ctparticles.select(ctpoffsets, ctpcounts).read(ctparts, xfer_props);
278  }
279  }
280 
281  private:
282 
283  std::vector<int> version;
284  std::vector<std::vector<double> > evts, parts, pinfo;
285  std::vector<std::vector<double> > ctevts, ctparts;
286  std::vector<std::string> wgtnames;
287 
288  inline Particle GetParticle(size_t i) const {
289  return Particle(parts[i][0], parts[i][1], parts[i][2], parts[i][3],
290  parts[i][4], parts[i][5], parts[i][6], parts[i][7],
291  parts[i][8], parts[i][9], parts[i][10],
292  parts[i][11], parts[i][12]);
293  }
294 
295  inline Particle GetCTParticle(size_t i) const {
296  return Particle(-1, -1, -1, -1, -1, -1, ctparts[i][0], ctparts[i][1],
297  ctparts[i][2], ctparts[i][3], -1, -1, -1);
298  }
299 
300 };
301 
302 //==========================================================================
303 
304 }// end of namespace LHEH5
305 
306 #endif // Pythia8_LHEH5v2_H
Definition: LHEH5.h:59
Standard includes.
Process info struct.
Definition: LHEH5v2.h:67
LHEFile class.
Definition: LHEH5v2.h:172
double UnitWeight() const
Definition: LHEH5v2.h:190
Event struct.
Definition: LHEH5v2.h:121
Particle struct.
Definition: LHEH5.h:65