PYTHIA  8.312
LHAHDF5v2.h
1 // LHAHDF5v2.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2024 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 // Authors: Christian Preuss, Stefan Hoeche December 2023.
6 
7 #ifndef Pythia8_LHAHDF5v2_H
8 #define Pythia8_LHAHDF5v2_H
9 
10 // Interface includes.
11 #include "Pythia8Plugins/LHEH5v2.h"
12 
13 // Generator includes.
14 #include "Pythia8/Pythia.h"
15 
16 namespace Pythia8 {
17 
18 //==========================================================================
19 
20 // HDF5 version 2 file reader.
21 // Converts to Pythia-internal events by acting as replacement
22 // Les Houches Event reader.
23 
24 class LHAupH5v2 : public Pythia8::LHAup {
25 
26  public:
27 
28  LHAupH5v2(HighFive::File* h5fileIn, size_t firstEventIn, size_t readSizeIn,
29  bool normalize) :
30  lhefPtr(new LHEH5::LHEFile()), readSizeSav(readSizeIn),
31  nReadSav(0), nTrialsSav(0) {
32 
33  // Read event-file header and events.
34  lhefPtr->ReadHeader(*h5fileIn);
35  lhefPtr->ReadEvents(*h5fileIn, firstEventIn, readSizeIn);
36 
37  // This reads the init information.
38  _weightnames = lhefPtr->WeightNames();
39  if (normalize) lhefPtr->Scale(readSizeIn/lhefPtr->SumTrials());
40  std::vector<double> info;
41  h5fileIn->getDataSet("init").read(info);
42  setBeamA(info[0],info[2],info[4],info[6]);
43  setBeamB(info[1],info[3],info[5],info[7]);
44  setStrategy(-4);
45  int numProcesses = info[9];
46  vector<int> procId(numProcesses);
47  vector<double> xSection(numProcesses);
48  vector<double> error(numProcesses);
49  vector<double> unitWeight(numProcesses);
50  for (int i = 0; i<numProcesses; ++i) {
51  LHEH5::ProcInfo pi(lhefPtr->GetProcInfo(i));
52  procId[i] = pi.pid;
53  xSection[i] = pi.xsec;
54  error[i] = pi.error;
55  unitWeight[i] = pi.unitwgt;
56  }
57  for (int np = 0; np<numProcesses; ++np) {
58  addProcess(procId[np], xSection[np], error[np], unitWeight[np]);
59  xSecSumSave += xSection[np];
60  xErrSumSave += pow2(error[np]);
61  }
62  }
63 
64  ~LHAupH5v2() {delete lhefPtr;}
65 
66  bool setInit() override {return true;}
67  bool setEvent(int idProc=0) override;
68  void forceStrategy(int strategyIn) {setStrategy(strategyIn);}
69  size_t nTrials() {return nTrialsSav;}
70  size_t nRead() {return nReadSav;}
71 
72 private:
73 
74  // HDF5 event file.
75  LHEH5::LHEFile *lhefPtr;
76 
77  // Info for reader.
78  size_t readSizeSav, nReadSav, nTrialsSav;
79 
80  // Multiweight vector.
81  vector<double> _eventweightvalues;
82  vector<string> _weightnames;
83 
84  // Particle production scales.
85  LHAscales scalesNow;
86 
87 };
88 
89 //--------------------------------------------------------------------------
90 
91 // Read an event.
92 
94 
95  // Equivalent of end of file.
96  if (nReadSav >= readSizeSav) return false;
97 
98  // Read event.
99  LHEH5::Event evt(lhefPtr->GetEvent(nReadSav));
100  if (evt[0].pz<0 && evt[1].pz>0) swap<LHEH5::Particle>(evt[0], evt[1]);
101 
102  setProcess(evt.pinfo.pid, evt.wgts[0], evt.mur, evt.aqed, evt.aqcd);
103  nupSave = evt.size();
104  idprupSave = evt.pinfo.pid;
105  xwgtupSave = evt.wgts[0];
106  scalupSave = evt.mur;
107  aqedupSave = evt.aqed;
108  aqcdupSave = evt.aqcd;
109  double scalein = -1.;
110 
111  // Communicate event weight to Info.
112  _eventweightvalues=evt.wgts;
113  // infoPtr->weights_compressed_names = &_weightnames;
114  infoPtr->weights_compressed = &_eventweightvalues;
115 
116  // Set particles.
117  for (unsigned int ip=0; ip<evt.size(); ++ip) {
118  const LHEH5::Particle& p = evt[ip];
119  if (ip < 2) addParticle(p.id, p.st, 0, 0,
120  p.cl1, p.cl2, p.px, p.py, p.pz, p.e, p.m,
121  p.lt, p.sp, scalein);
122  else addParticle(p.id, p.st, p.mo1, p.mo2,
123  p.cl1, p.cl2, p.px, p.py, p.pz, p.e, p.m,
124  p.lt, p.sp, scalein);
125  }
126 
127  // Scale setting
128  scalesNow.clear();
129  scalesNow.muf = evt.muf;
130  scalesNow.mur = evt.mur;
131  scalesNow.mups = evt.muq;
132  infoPtr->scales = &scalesNow;
133  infoPtr->setEventAttribute("npLO",std::to_string(evt.pinfo.nplo));
134  infoPtr->setEventAttribute("npNLO",std::to_string(evt.pinfo.npnlo));
135 
136  // Update counters.
137  ++nReadSav;
138  nTrialsSav += evt.trials;
139  return true;
140 }
141 
142 //==========================================================================
143 
144 } // end namespace Pythia8
145 
146 #endif // Pythia8_LHAHDF5v2_H
void setBeamA(int idIn, double eIn, int pdfGroupIn=0, int pdfSetIn=0)
Input beam info.
Definition: LesHouches.h:222
void addProcess(int idProcIn, double xSecIn=1., double xErrIn=0., double xMaxIn=1.)
Input process info.
Definition: LesHouches.h:233
constexpr double pow2(const double &x)
Powers of small integers - for balance speed/code clarity.
Definition: PythiaStdlib.h:182
Definition: LHEH5.h:59
void setStrategy(int strategyIn)
Input process weight strategy.
Definition: LesHouches.h:230
vector< double > * weights_compressed
The weights associated with this event, as given by the LHAweights tags.
Definition: Info.h:392
double mups
Definition: LHEF3.h:308
Collect different scales relevant for an event.
Definition: LHEF3.h:281
Process info struct.
Definition: LHEH5v2.h:67
void setEventAttribute(string key, string value, bool doOverwrite=true)
Externally set event tag auxiliary information.
Definition: Info.h:431
Definition: LesHouches.h:78
LHAscales * scales
Contents of the LHAscales tag.
Definition: Info.h:395
void addParticle(LHAParticle particleIn)
Input particle info, one particle at the time.
Definition: LesHouches.h:252
void setProcess(int idProcIn=0, double weightIn=1., double scaleIn=0., double alphaQEDIn=0.0073, double alphaQCDIn=0.12)
Input info on the selected process.
Definition: LesHouches.h:243
Definition: LHAHDF5v2.h:24
LHEFile class.
Definition: LHEH5v2.h:172
double muf
The factorization scale used for this event.
Definition: LHEF3.h:301
void clear()
Function to reset this object.
Definition: LHEF3.h:294
int nupSave
Event properties from LHEF files, for repeated use.
Definition: LesHouches.h:298
bool setInit() override
Definition: LHAHDF5v2.h:66
double mur
The renormalization scale used for this event.
Definition: LHEF3.h:304
Header for classes to set beam momentum and interaction vertex spread.
Definition: Analysis.h:20
Info * infoPtr
Pointer to various information on the generation.
Definition: LesHouches.h:218
bool setEvent(int idProc=0) override
Read an event.
Definition: LHAHDF5v2.h:93
Event struct.
Definition: LHEH5v2.h:121
LHAupH5v2(HighFive::File *h5fileIn, size_t firstEventIn, size_t readSizeIn, bool normalize)
Definition: LHAHDF5v2.h:28
Particle struct.
Definition: LHEH5.h:65