PYTHIA  8.312
LHAHDF5.h
1 // LHAHDF5.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 // Author: Christian Preuss, January 2021.
6 
7 #ifndef Pythia8_LHAHDF5_H
8 #define Pythia8_LHAHDF5_H
9 
10 // Interface includes.
11 #include "Pythia8Plugins/LHEH5.h"
12 
13 // Generator includes.
14 #include "Pythia8/Pythia.h"
15 
16 namespace Pythia8 {
17 
18 //==========================================================================
19 
20 // HDF5 file reader. Converts to Pythia-internal events by acting as
21 // replacement Les Houches Event reader.
22 
23 class LHAupH5 : public Pythia8::LHAup {
24 
25  public:
26 
27  LHAupH5(HighFive::File* h5fileIn, size_t firstEventIn, size_t readSizeIn,
28  string versionIn = "") :
29  nTrialsSav(0), nRead(0), isOldFormat(h5fileIn->exist("/index")),
30  particleSav(h5fileIn->getGroup("particle")),
31  eventSav(h5fileIn->getGroup("event")),
32  initSav(h5fileIn->getGroup("init")),
33  procInfoSav(h5fileIn->getGroup("procInfo")),
34  readSizeSav(readSizeIn), firstEventSav(firstEventIn), fileSav(h5fileIn) {
35 
36  // Check if version number exists in event file.
37  valid = false;
38  if (h5fileIn->exist("/init/version")) {
39  DataSet version = initSav.getDataSet("version");
40  vector<int> versionNo;
41  version.read(versionNo);
42  versionSav = to_string(versionNo[0]) + "." + to_string(versionNo[1])
43  + "." + to_string(versionNo[2]);
44  isOldFormat = (versionSav=="0.1.0");
45  } else if (isOldFormat) {
46  // Based on the existence of the index group, we can guess this is 0.1.0.
47  versionSav = "0.1.0";
48  } else {
49  // Have to rely on input if we could not guess the version.
50  versionSav = versionIn;
51  if (versionSav=="") {
52  cout << " LHAupH5 error - cannot determine version number.\n";
53  return;
54  }
55  }
56  cout << " LHAupH5 version " << versionSav << ".\n";
57 
58  // Check if version supports multiweights. (Starting from 1.0.0.)
59  hasMultiWts = !(versionSav=="0.1.0" || versionSav=="0.2.0");
60  cout << " LHAupH5 file format "
61  << (hasMultiWts?"supports":"does not support")
62  << " multi-weights." << endl;
63 
64  hid_t dspace;
65  if( !isOldFormat ) {
66  DataSet npLO = procInfoSav.getDataSet("npLO");
67  DataSet npNLO = procInfoSav.getDataSet("npNLO");
68  npLO.read(npLOSav);
69  npNLO.read(npNLOSav);
70  dspace = H5Dget_space(fileSav->getDataSet("event/start").getId());
71  } else {
72  indexSav = fileSav->getGroup("index");
73  dspace = H5Dget_space(fileSav->getDataSet("index/start").getId());
74  }
75  // Check if size of read is compatible.
76  if (readSizeIn > H5Sget_simple_extent_npoints(dspace) ) {
77  cout << "H5 size request incompatible with file.\n";
78  return;
79  }
80 
81  // Check for multiweights.
82  if ( fileSav->exist("event/weight") && hasMultiWts ) {
83  DataSet weights = eventSav.getDataSet("weight");
84  auto attr_keys = weights.listAttributeNames();
85  Attribute a = weights.getAttribute(attr_keys[0]);
86  a.read(weightsNames);
87  }
88 
89  // This reads and holds the information of readSize events,
90  // starting from firstEvent.
91  if (!isOldFormat) {
92  lheEvts2Sav = LHEH5::readEvents2(particleSav, eventSav,
93  firstEventSav, readSizeSav, npLOSav, npNLOSav, hasMultiWts);
94 
95  } else lheEvtsSav = LHEH5::readEvents(indexSav, particleSav, eventSav,
96  firstEventSav, readSizeSav);
97  valid = true;
98 
99  }
100 
101  // Read and set the info from init and procInfo.
102  bool setInit() override;
103  bool setEvent(int idProc=0) override;
104  void forceStrategy(int strategyIn) {setStrategy(strategyIn);}
105  size_t nTrials() {return nTrialsSav;}
106 
107  private:
108 
109  // HDF5 file.
110  HighFive::File* fileSav{};
111 
112  // HDF5 Groups.
113  HighFive::Group indexSav, particleSav, eventSav, initSav, procInfoSav;
114 
115  // Events from HDF5 file.
116  LHEH5::Events lheEvtsSav;
117  LHEH5::Events2 lheEvts2Sav;
118 
119  // Info for reader.
120  size_t readSizeSav, firstEventSav, nTrialsSav;
121  int npLOSav, npNLOSav;
122  bool valid, isOldFormat, hasMultiWts;
123  string versionSav;
124 
125  // Additional parameters.
126  int nRead;
127 
128  // Multiweight vector. Reset each event.
129  vector<double> weightsSav;
130  vector<string> weightsNames;
131 
132  // Particle production scales.
133  LHAscales scalesNow;
134 
135 };
136 
137 //--------------------------------------------------------------------------
138 
139 // HDF5 file reader. Converts to Pythia-internal events by acting as
140 // replacement Les Houches Event reader.
141 
143  int beamA, beamB;
144  double energyA, energyB;
145  int PDFgroupA, PDFgroupB;
146  int PDFsetA, PDFsetB;
147 
148  if (!valid) return false;
149  initSav.getDataSet("beamA").read(beamA);
150  initSav.getDataSet("energyA").read(energyA);
151  initSav.getDataSet("PDFgroupA").read(PDFgroupA);
152  initSav.getDataSet("PDFsetA").read(PDFsetA);
153 
154  initSav.getDataSet("beamB").read(beamB);
155  initSav.getDataSet("energyB").read(energyB);
156  initSav.getDataSet("PDFgroupB").read(PDFgroupB);
157  initSav.getDataSet("PDFsetB").read(PDFsetB);
158 
159  setBeamA(beamA, energyA, PDFgroupA, PDFsetA);
160  setBeamB(beamB, energyB, PDFgroupB, PDFsetB);
161 
162  int weightingStrategy = 3;
163  initSav.getDataSet("weightingStrategy").read(weightingStrategy);
164  setStrategy(weightingStrategy);
165 
166  int nProcesses = 1;
167  initSav.getDataSet("numProcesses").read(nProcesses);
168 
169  vector<int> procId;
170  vector<double> xSection, error, unitWeight;
171  procInfoSav.getDataSet("procId").read(procId);
172  procInfoSav.getDataSet("xSection").read(xSection);
173  procInfoSav.getDataSet("error").read(error);
174  procInfoSav.getDataSet("unitWeight").read(unitWeight);
175  infoPtr->sigmaLHEFSave.resize(0);
176  for (int iProc(0); iProc<nProcesses; ++iProc) {
177  addProcess(procId[iProc], xSection[iProc], error[iProc],
178  unitWeight[iProc]);
179  infoPtr->sigmaLHEFSave.push_back(xSection[iProc]);
180  }
181  return true;
182 
183 }
184 
185 //--------------------------------------------------------------------------
186 
187 // Read an event.
188 
189 bool LHAupH5::setEvent(int idProc) {
190 
191  // Equivalent of end of file.
192  if (!valid) return false;
193  if (nRead >= readSizeSav) return false;
194  LHEH5::EventHeader evtHeader = !isOldFormat ?
195  lheEvts2Sav.mkEventHeader(nRead) : lheEvtsSav.mkEventHeader(nRead);
196  weightsSav = evtHeader.weights;
197  nTrialsSav += evtHeader.trials;
198  // Skip zero-weight events, but add trials.
199  while (weightsSav[0] == 0. && nRead < readSizeSav - 1) {
200  ++nRead;
201  evtHeader = !isOldFormat ?
202  lheEvts2Sav.mkEventHeader(nRead) : lheEvtsSav.mkEventHeader(nRead);
203  weightsSav = evtHeader.weights;
204  nTrialsSav += evtHeader.trials;
205  }
206  // Communicate event weight to Info.
207  infoPtr->weightContainerPtr->setWeightNominal( weightsSav[0] );
208  infoPtr->weightContainerPtr->weightsLHEF.bookVectors(
209  weightsSav, weightsNames );
210 
211  xwgtupSave = weightsSav[0];
212  idprupSave = evtHeader.pid;
213  scalupSave = evtHeader.scale;
214  aqedupSave = evtHeader.aqed;
215  aqcdupSave = evtHeader.aqcd;
216  setProcess(idprupSave, xwgtupSave, scalupSave, aqedupSave, aqcdupSave);
217  double scalein = -1.;
218  vector<LHEH5::Particle> particles;
219  particles = !isOldFormat ? lheEvts2Sav.mkEvent(nRead)
220  : lheEvtsSav.mkEvent(nRead);
221 
222  // Set particles.
223  int nPtcls = 0;
224  for (int iPtcl(0); iPtcl<int(particles.size()); ++iPtcl) {
225  LHEH5::Particle ptcl = particles.at(iPtcl);
226  if (ptcl.id == 0) continue;
227  nPtcls++;
228  if (iPtcl < 2) addParticle(ptcl.id, ptcl.status, 0, 0,
229  ptcl.color1, ptcl.color2, ptcl.px, ptcl.py, ptcl.pz, ptcl.e, ptcl.m,
230  ptcl.lifetime, ptcl.spin, scalein);
231  else addParticle(ptcl.id, ptcl.status, ptcl.mother1, ptcl.mother2,
232  ptcl.color1, ptcl.color2, ptcl.px, ptcl.py, ptcl.pz, ptcl.e, ptcl.m,
233  ptcl.lifetime, ptcl.spin, scalein);
234  }
235  nupSave = nPtcls;
236 
237  // Scale setting
238  scalesNow.clear();
239  scalesNow.muf = evtHeader.fscale;
240  scalesNow.mur = evtHeader.rscale;
241  scalesNow.mups = evtHeader.scale;
242  infoPtr->scales = &scalesNow;
243  ++nRead;
244  return true;
245 
246 }
247 
248 //==========================================================================
249 
250 } // end namespace Pythia8
251 
252 #endif // Pythia8_LHAHDF5_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
EventHeader mkEventHeader(int ievent) const
Make an event header, given an index.
Definition: LHEH5.h:423
Old event struct.
Definition: LHEH5.h:127
void bookVectors(vector< double > weights, vector< string > names) override
Store the current event information.
Definition: Weights.cc:360
Definition: LHAHDF5.h:23
void setStrategy(int strategyIn)
Input process weight strategy.
Definition: LesHouches.h:230
WeightsLHEF weightsLHEF
First example of a weight subcategory.
Definition: Weights.h:408
EventHeader mkEventHeader(int ievent) const
Make an event header, given an index.
Definition: LHEH5.h:200
std::vector< double > weights
this is all LHAu-::setProcess
Definition: LHEH5.h:96
New events struct.
Definition: LHEH5.h:351
double mups
Definition: LHEF3.h:308
Collect different scales relevant for an event.
Definition: LHEF3.h:281
std::vector< Particle > mkEvent(size_t ievent) const
Make an event, given an index.
Definition: LHEH5.h:180
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
bool setEvent(int idProc=0) override
Read an event.
Definition: LHAHDF5.h:189
LHAupH5(HighFive::File *h5fileIn, size_t firstEventIn, size_t readSizeIn, string versionIn="")
Definition: LHAHDF5.h:27
void setWeightNominal(double weightNow)
The WeightContainer class.
Definition: Weights.cc:852
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
std::vector< Particle > mkEvent(size_t ievent) const
Make an event, given an index.
Definition: LHEH5.h:404
double mur
The renormalization scale used for this event.
Definition: LHEF3.h:304
size_t trials
double weight;
Definition: LHEH5.h:98
bool setInit() override
Read and set the info from init and procInfo.
Definition: LHAHDF5.h:142
Header for classes to set beam momentum and interaction vertex spread.
Definition: Analysis.h:20
int pid
corr to NUP
Definition: LHEH5.h:95
Info * infoPtr
Pointer to various information on the generation.
Definition: LesHouches.h:218
Event header struct.
Definition: LHEH5.h:92
vector< double > sigmaLHEFSave
Save process information as read from init block of LHEF.
Definition: Info.h:370
Particle struct.
Definition: LHEH5.h:65