Skip to content

Commit 2f36720

Browse files
cholmcccholmcc
committed
Full header read-in and external program
This change does two things: **Full header** _All_ information available in the HepMC event header is propagated to the MC event header information map. This includes - Heavy-ion "geometry" parameters (b,Ncoll,Npart,...) - Cross-section(s) - Weight(s) - PDF information - and other attributes defined This is so that we can build a full HepMC event structure later - for example to pass to Rivet analyses ** External program ** The functionality of the generatator is expanded so that it may spawn an event generator program, say `eg`. - The generator opens a FIFO - The generator then executes the program `eg` in the background - The `eg` program is assumed to write HepMC event records on standard output, which is then redirected to the FIFO - The generator reads events from the FIFO For this to work, a number of conditions _must_ be met by the `eg` program: - It _must_ write events in the HepMC event format - It _must_ write the HepMC event records to standard output - It _cannot_ write anything else but the HepMC event record to standard output - It _must_ accept the command line option `-n NEVENTS` to set the number of events to generate. If a particular `eg` program does not meet these requirements, then a simple shell script can be defined to wrap the `eg` appropriately. For example, the CRMC program `crmc` _can_ write HepMC events to standard output, but it will also dump other stuff there. Thus, we can provide the script #!/bin/sh crmc $@ -o hepmc3 -f /dev/stdout | \ sed -n 's/^\(HepMC::\|[EAUWVP] \)/\1/p' which simply filters the output of `crmc`. Another EG program may not accept the `-n EVENTS` command line option, but rather has the command line option `--nevents`, so then we would do something like #!/bin/sh cmdline="eg-program -o /dev/stdout " while test $# -gt 0 ; do case x$1 in x-n) cmdline="$cmdline -n $2"; shift ;; *) cmdline="$cmdline $1" ;; esac shift done $cmdline The command line to run is specified as --configKeyValues "HepMC.progCmd=<program and options>" and can include not only the program name but also other options to the program. For example --configKeyValues "HepMC.progCmd=crmc -m 5 -i 20800820 -I 20800820" for Pb-Pb collisions with Hijing. With this change, we can use _any_ event generator which is capable to write out its event records in the HepMC format.
1 parent dadc5c6 commit 2f36720

3 files changed

Lines changed: 198 additions & 2 deletions

File tree

‎Generators/include/Generators/GeneratorHepMC.h‎

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -62,8 +62,9 @@ class GeneratorHepMC : public Generator
6262
/** setters **/
6363
void setVersion(Int_t val) { mVersion = val; };
6464
void setFileName(std::string val) { mFileName = val; };
65+
void setProgCmd(std::string val) { mProgCmd = val; };
6566
void setEventsToSkip(uint64_t val) { mEventsToSkip = val; };
66-
67+
void setNEvents(unsigned int val) { mNEvents = val; }
6768
protected:
6869
/** copy constructor **/
6970
GeneratorHepMC(const GeneratorHepMC&);
@@ -77,11 +78,16 @@ class GeneratorHepMC : public Generator
7778
const HepMC3::FourVector getBoostedVector(const HepMC3::FourVector& vector, Double_t boost);
7879
#endif
7980

81+
/** methods that can be overridded **/
82+
void updateHeader(o2::dataformats::MCEventHeader* eventHeader) override;
83+
8084
/** HepMC interface **/
8185
std::ifstream mStream; //!
8286
std::string mFileName;
87+
std::string mProgCmd;
8388
Int_t mVersion;
8489
uint64_t mEventsToSkip;
90+
unsigned int mNEvents;
8591
#ifdef GENERATORS_WITH_HEPMC3_DEPRECATED
8692
HepMC::Reader* mReader; //!
8793
HepMC::GenEvent* mEvent; //!

‎Generators/include/Generators/GeneratorHepMCParam.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -31,6 +31,7 @@ namespace eventgen
3131

3232
struct GeneratorHepMCParam : public o2::conf::ConfigurableParamHelper<GeneratorHepMCParam> {
3333
std::string fileName = "";
34+
std::string progCmd = "";//Program command line to spawn, must write HepMC on stdout
3435
int version = 2;
3536
uint64_t eventsToSkip = 0;
3637
O2ParamDef(GeneratorHepMCParam, "HepMC");

‎Generators/src/GeneratorHepMC.cxx‎

Lines changed: 190 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -14,15 +14,22 @@
1414
#include "SimulationDataFormat/MCUtils.h"
1515
#include "Generators/GeneratorHepMC.h"
1616
#include "Generators/GeneratorHepMCParam.h"
17+
#include "SimulationDataFormat/MCEventHeader.h"
1718
#include "HepMC3/ReaderAscii.h"
1819
#include "HepMC3/ReaderAsciiHepMC2.h"
1920
#include "HepMC3/GenEvent.h"
2021
#include "HepMC3/GenParticle.h"
2122
#include "HepMC3/GenVertex.h"
2223
#include "HepMC3/FourVector.h"
24+
#include "HepMC3/Version.h"
2325
#include "TParticle.h"
2426
#include "TSystem.h"
2527

28+
#include <cstdlib>
29+
#include <sys/types.h> // POSIX only
30+
#include <sys/stat.h> // POISX only
31+
#include <cstdio>
32+
2633
#include <fairlogger/Logger.h>
2734
#include "FairPrimaryGenerator.h"
2835
#include <cmath>
@@ -85,6 +92,7 @@ Bool_t GeneratorHepMC::generateEvent()
8592
mEvent->clear();
8693
mReader->read_event(*mEvent);
8794
if (mReader->failed()) {
95+
LOG(error) << "Failed to read one event from input";
8896
return kFALSE;
8997
}
9098
/** set units to desired output **/
@@ -143,6 +151,133 @@ Bool_t GeneratorHepMC::importParticles()
143151
return kTRUE;
144152
}
145153

154+
namespace
155+
{
156+
void putAttributeInfo(o2::dataformats::MCEventHeader* eventHeader,
157+
const std::string& name,
158+
const std::shared_ptr<HepMC3::Attribute>& a)
159+
{
160+
if (auto* p = dynamic_cast<HepMC3::IntAttribute*>(a.get()))
161+
eventHeader->putInfo<int>(name,p->value());
162+
if (auto* p = dynamic_cast<HepMC3::LongAttribute*>(a.get()))
163+
eventHeader->putInfo<int>(name,p->value());
164+
if (auto* p = dynamic_cast<HepMC3::FloatAttribute*>(a.get()))
165+
eventHeader->putInfo<float>(name,p->value());
166+
if (auto* p = dynamic_cast<HepMC3::DoubleAttribute*>(a.get()))
167+
eventHeader->putInfo<float>(name,p->value());
168+
if (auto* p = dynamic_cast<HepMC3::StringAttribute*>(a.get()))
169+
eventHeader->putInfo<std::string>(name,p->value());
170+
if (auto* p = dynamic_cast<HepMC3::CharAttribute*>(a.get()))
171+
eventHeader->putInfo<char>(name,p->value());
172+
if (auto* p = dynamic_cast<HepMC3::LongLongAttribute*>(a.get()))
173+
eventHeader->putInfo<int>(name,p->value());
174+
if (auto* p = dynamic_cast<HepMC3::LongDoubleAttribute*>(a.get()))
175+
eventHeader->putInfo<float>(name,p->value());
176+
if (auto* p = dynamic_cast<HepMC3::UIntAttribute*>(a.get()))
177+
eventHeader->putInfo<int>(name,p->value());
178+
if (auto* p = dynamic_cast<HepMC3::ULongAttribute*>(a.get()))
179+
eventHeader->putInfo<int>(name,p->value());
180+
if (auto* p = dynamic_cast<HepMC3::ULongLongAttribute*>(a.get()))
181+
eventHeader->putInfo<int>(name,p->value());
182+
if (auto* p = dynamic_cast<HepMC3::BoolAttribute*>(a.get()))
183+
eventHeader->putInfo<bool>(name,p->value());
184+
}
185+
}
186+
187+
188+
/*****************************************************************/
189+
190+
void GeneratorHepMC::updateHeader(o2::dataformats::MCEventHeader* eventHeader)
191+
{
192+
/** update header **/
193+
using Key=o2::dataformats::MCInfoKeys;
194+
195+
eventHeader->putInfo<std::string>(Key::generator, "hepmc");
196+
eventHeader->putInfo<int>(Key::generatorVersion, HEPMC3_VERSION_CODE);
197+
198+
auto xSection = mEvent->cross_section();
199+
auto pdfInfo = mEvent->pdf_info();
200+
auto hiInfo = mEvent->heavy_ion();
201+
202+
// Set default cross-section
203+
if (xSection) {
204+
eventHeader->putInfo<float>(Key::xSection, xSection->xsec());
205+
eventHeader->putInfo<float>(Key::xSectionError, xSection->xsec_err());
206+
eventHeader->putInfo<int>(Key::acceptedEvents,
207+
xSection->get_accepted_events());
208+
eventHeader->putInfo<int>(Key::attemptedEvents,
209+
xSection->get_attempted_events());
210+
}
211+
212+
// Set weights and cross sections
213+
size_t iw = 0;
214+
for (auto w : mEvent->weights()) {
215+
std::string post = (iw > 0 ? "_"+std::to_string(iw) : "");
216+
eventHeader->putInfo<float>(Key::weight+post, w);
217+
if (xSection) {
218+
eventHeader->putInfo<float>(Key::xSection, xSection->xsec(iw));
219+
eventHeader->putInfo<float>(Key::xSectionError,xSection->xsec_err(iw));
220+
}
221+
iw++;
222+
}
223+
224+
// Set the PDF information
225+
if (pdfInfo) {
226+
eventHeader->putInfo<int>(Key::pdfParton1Id, pdfInfo->parton_id[0]);
227+
eventHeader->putInfo<int>(Key::pdfParton2Id, pdfInfo->parton_id[1]);
228+
eventHeader->putInfo<float>(Key::pdfX1, pdfInfo->x[0]);
229+
eventHeader->putInfo<float>(Key::pdfX2, pdfInfo->x[1]);
230+
eventHeader->putInfo<float>(Key::pdfScale, pdfInfo->scale);
231+
eventHeader->putInfo<float>(Key::pdfXF1, pdfInfo->xf[0]);
232+
eventHeader->putInfo<float>(Key::pdfXF2, pdfInfo->xf[1]);
233+
eventHeader->putInfo<float>(Key::pdfCode1, pdfInfo->pdf_id[0]);
234+
eventHeader->putInfo<float>(Key::pdfCode2, pdfInfo->pdf_id[1]);
235+
}
236+
237+
// Set heavy-ion information
238+
if (hiInfo) {
239+
eventHeader->putInfo<int>(Key::impactParameter ,
240+
hiInfo->impact_parameter);
241+
eventHeader->putInfo<int>(Key::nPart ,
242+
hiInfo->Npart_proj+hiInfo->Npart_targ);
243+
eventHeader->putInfo<int>(Key::nPartProjectile ,hiInfo->Npart_proj);
244+
eventHeader->putInfo<int>(Key::nPartTarget ,hiInfo->Npart_targ);
245+
eventHeader->putInfo<int>(Key::nColl ,hiInfo->Ncoll);
246+
eventHeader->putInfo<int>(Key::nCollHard ,hiInfo->Ncoll_hard);
247+
eventHeader->putInfo<int>(Key::nCollNNWounded ,
248+
hiInfo->N_Nwounded_collisions);
249+
eventHeader->putInfo<int>(Key::nCollNWoundedN ,
250+
hiInfo->Nwounded_N_collisions);
251+
eventHeader->putInfo<int>(Key::nCollNWoundedNwounded ,
252+
hiInfo->Nwounded_Nwounded_collisions);
253+
eventHeader->putInfo<int>(Key::planeAngle ,
254+
hiInfo->event_plane_angle);
255+
eventHeader->putInfo<int>(Key::sigmaInelNN ,
256+
hiInfo->sigma_inel_NN);
257+
eventHeader->putInfo<int>(Key::centrality ,hiInfo->centrality);
258+
eventHeader->putInfo<int>(Key::nSpecProjectileProton ,hiInfo->Nspec_proj_p);
259+
eventHeader->putInfo<int>(Key::nSpecProjectileNeutron,hiInfo->Nspec_proj_n);
260+
eventHeader->putInfo<int>(Key::nSpecTargetProton ,hiInfo->Nspec_targ_p);
261+
eventHeader->putInfo<int>(Key::nSpecTargetNeutron ,hiInfo->Nspec_targ_n);
262+
}
263+
264+
for (auto na : mEvent->attributes()) {
265+
std::string name = na.first;
266+
if (name == "GenPdfInfo" ||
267+
name == "GenCrossSection" ||
268+
name == "GenHeavyIon") continue;
269+
270+
for (auto ia : na.second) {
271+
int no = ia.first;
272+
auto at = ia.second;
273+
std::string post = (no == 0 ? "" : std::to_string(no));
274+
275+
putAttributeInfo(eventHeader, name+post, at);
276+
}
277+
}
278+
}
279+
280+
146281
/*****************************************************************/
147282

148283
Bool_t GeneratorHepMC::Init()
@@ -152,14 +287,68 @@ Bool_t GeneratorHepMC::Init()
152287
/** init base class **/
153288
Generator::Init();
154289

155-
/** open file **/
156290
std::string filename = gSystem->ExpandPathName(mFileName.c_str());
291+
292+
// If a EG command line is given, then we make a fifo on a temporary
293+
// file, and directs the EG to write to that fifo. We will then set
294+
// up the HepMC3 reader to read from that fifo.
295+
//
296+
// o2-sim -g hepmc --configKeyValues "HepMC.progCmd=<cmd>" ...
297+
//
298+
// where <cmd> is the command line to run an event generator. The
299+
// event generator should output HepMC event records to standard
300+
// output. Nothing else, but the HepMC event record may be output
301+
// to standard output. If the EG has other output to standard
302+
// output, then a filter can be set-up. For example
303+
//
304+
// crmc -n 3 -o hepmc3 -c /optsw/inst/etc/crmc.param -f /dev/stdout \
305+
// | sed -n 's/^\(HepMC::\|[EAUWVP] \)/\1/p'
306+
//
307+
// What's more, the event generator program must accept the command
308+
// line argument `-n NEVENTS` to set the number of events to
309+
// produce.
310+
//
311+
// All of this can conviniently be achieved via a wrapper script
312+
// around the actual EG program.
313+
if (not mProgCmd.empty()) {
314+
// Set filename to be a temporary name
315+
// Should perhaps use
316+
//
317+
// TString base("xxxxxx");
318+
// auto fp = gSystem->TempFileName(base);
319+
// fclose(fp);
320+
//
321+
filename = std::tmpnam(nullptr);
322+
323+
// Make a fifo
324+
int ret = mkfifo(filename.c_str(), 0600);
325+
if (ret != 0) {
326+
LOG(fatal) << "Failed to make fifo \"" << filename << "\"";
327+
return false;
328+
}
329+
330+
// Build command line, rediret stdout to our fifo and put
331+
// in the background.
332+
std::string cmd =
333+
mProgCmd +
334+
" -n " + std::to_string(mNEvents) +
335+
" > " + filename + " &";
336+
LOG(info) << "EG command line is \"" << cmd << "\"";
337+
338+
ret = std::system(cmd.c_str());
339+
if (ret != 0) {
340+
LOG(fatal) << "Failed to spawn \"" << cmd << "\"";
341+
return false;
342+
}
343+
}
344+
/** open file **/
157345
mStream.open(filename);
158346
if (!mStream.is_open()) {
159347
LOG(fatal) << "Cannot open input file: " << filename << std::endl;
160348
return kFALSE;
161349
}
162350

351+
LOG(info) << "Set up reader to read from \"" << filename << "\"" << std::endl;
163352
/** create reader according to HepMC version **/
164353
switch (mVersion) {
165354
case 2:

0 commit comments

Comments
 (0)