Skip to content

Commit 0b8c28d

Browse files
jackal1-66jackal1-66
authored
Herwig7 external generator (#2163)
Implementation of HERWIG7 external generator to simplify usage
1 parent a6fa4ed commit 0b8c28d

File tree

4 files changed

+410
-0
lines changed

4 files changed

+410
-0
lines changed

‎MC/config/examples/herwig/LHC.in‎

Lines changed: 36 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,36 @@
1+
# -*- ThePEG-repository -*-
2+
3+
################################################################################
4+
# This file contains our best tune to UE data from ATLAS at 7 TeV. More recent
5+
# tunes and tunes for other centre-of-mass energies as well as more usage
6+
# instructions can be obtained from this Herwig wiki page:
7+
# http://projects.hepforge.org/herwig/trac/wiki/MB_UE_tunes
8+
# The model for soft interactions and diffractions is explained in
9+
# [S. Gieseke, P. Kirchgaesser, F. Loshaj, arXiv:1612.04701]
10+
################################################################################
11+
12+
read snippets/PPCollider.in
13+
14+
##################################################
15+
# Technical parameters for this run
16+
##################################################
17+
cd /Herwig/Generators
18+
##################################################
19+
# LHC physics parameters (override defaults here)
20+
##################################################
21+
set EventGenerator:EventHandler:LuminosityFunction:Energy 13600.0
22+
# Minimum Bias
23+
read snippets/MB.in
24+
25+
# Recommended set of parameters for MB/UE runs
26+
set /Herwig/Hadronization/ColourReconnector:ReconnectionProbability 0.5
27+
set /Herwig/UnderlyingEvent/MPIHandler:pTmin0 3.502683
28+
set /Herwig/UnderlyingEvent/MPIHandler:InvRadius 1.402055
29+
set /Herwig/UnderlyingEvent/MPIHandler:Power 0.416852
30+
set /Herwig/Partons/RemnantDecayer:ladderPower -0.08
31+
set /Herwig/Partons/RemnantDecayer:ladderNorm 0.95
32+
33+
##################################################
34+
# Save run for later usage with 'Herwig run'
35+
##################################################
36+
saverun LHC EventGenerator

‎MC/config/examples/herwig/generator_Herwig.C‎

Lines changed: 333 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,333 @@
1+
#define SKIP_HEPMC_CONVERSION 1
2+
#define HAVE_HEPMC3 1
3+
4+
// O2DPG and ROOT includes
5+
#include "FairGenerator.h"
6+
#include "FairPrimaryGenerator.h"
7+
#include "fairlogger/Logger.h"
8+
#include "TParticle.h"
9+
#include "TLorentzVector.h"
10+
#include <iostream>
11+
#include <vector>
12+
#include <string>
13+
14+
// ThePEG
15+
#include "ThePEG/Repository/EventGenerator.h"
16+
#include "ThePEG/EventRecord/Event.h"
17+
#include "ThePEG/EventRecord/Particle.h"
18+
#include "ThePEG/EventRecord/Step.h"
19+
#include "ThePEG/Config/ThePEG.h"
20+
#include "ThePEG/PDT/ParticleData.h"
21+
#include "ThePEG/Repository/Repository.h"
22+
23+
// Herwig
24+
#include "Herwig/API/HerwigAPI.h"
25+
#include "Herwig/API/HerwigUI.h"
26+
27+
// Subclass of HerwigUI to provide minimal implementation of the abstract class
28+
class SimpleHerwigUI : public Herwig::HerwigUI
29+
{
30+
public:
31+
SimpleHerwigUI(const std::string &inFile,
32+
Herwig::RunMode::Mode mode = Herwig::RunMode::READ, int seed = 0)
33+
: m_inFile(inFile), m_mode(mode),
34+
m_in(inFile), m_out(std::cout), m_err(std::cerr), mSeed(seed)
35+
{
36+
if (!m_in)
37+
{
38+
LOG(fatal) << "Cannot open Herwig input file: " << inFile;
39+
exit(1);
40+
}
41+
std::string hDir = std::getenv("HERWIG_ROOT");
42+
if (!hDir.empty())
43+
{
44+
Dirs.push_back(hDir + "/share/Herwig");
45+
}
46+
}
47+
48+
Herwig::RunMode::Mode runMode() const override { return m_mode; }
49+
50+
std::string repository() const override {
51+
std::string rpo = std::getenv("HERWIG_ROOT");
52+
rpo.append("/share/Herwig/HerwigDefaults.rpo");
53+
return rpo;
54+
}
55+
std::string inputfile() const override { return m_inFile; }
56+
std::string setupfile() const override { return ""; }
57+
58+
bool resume() const override { return false; }
59+
bool tics() const override { return false; }
60+
std::string tag() const override { return ""; }
61+
std::string integrationList() const override { return ""; }
62+
63+
const std::vector<std::string> &prependReadDirectories() const override
64+
{
65+
return Dirs;
66+
}
67+
const std::vector<std::string> &appendReadDirectories() const override
68+
{
69+
static std::vector<std::string> empty;
70+
return empty;
71+
}
72+
73+
long N() const override { return 1; } // number of events
74+
int seed() const override { return mSeed; }
75+
int jobs() const override { return 1; }
76+
unsigned int jobSize() const override { return 1; }
77+
unsigned int maxJobs() const override { return 1; }
78+
79+
void quitWithHelp() const override { std::exit(1); }
80+
void quit() const override { std::exit(1); }
81+
82+
std::ostream &outStream() const override { return m_out; }
83+
std::ostream &errStream() const override { return m_err; }
84+
std::istream &inStream() const override { return m_in; }
85+
86+
private:
87+
std::string m_inFile;
88+
Herwig::RunMode::Mode m_mode;
89+
mutable std::ifstream m_in;
90+
std::ostream &m_out;
91+
std::ostream &m_err;
92+
std::vector<std::string> Dirs;
93+
int mSeed = 0;
94+
};
95+
96+
namespace o2
97+
{
98+
namespace eventgen
99+
{
100+
101+
/// HERWIG7 event generator using ThePEG interface
102+
/// Author: Marco Giacalone (marco.giacalone@cern.ch)
103+
/// Based on the O2DPG external generator configurations
104+
class GeneratorHerwig : public Generator
105+
{
106+
public:
107+
/// Default constructor
108+
GeneratorHerwig(const std::string& configFile = "LHC.in", int seed = -1)
109+
: mConfigFile(configFile)
110+
, mEventGenerator(nullptr)
111+
{
112+
LOG(info) << "HERWIG7 Generator construction";
113+
LOG(info) << "Config file: " << mConfigFile;
114+
std::string extension = mConfigFile.substr(mConfigFile.find_last_of("."));
115+
if( extension == ".in" ) {
116+
mIsInputFile = true;
117+
LOG(info) << "Using input file for configuration";
118+
} else if(extension == ".run") {
119+
mIsInputFile = false;
120+
LOG(info) << "Using run file for configuration";
121+
} else {
122+
LOG(fatal) << "No file extension found in config file: " << mConfigFile;
123+
exit(1);
124+
}
125+
if (seed < 0)
126+
{
127+
auto &conf = o2::conf::SimConfig::Instance();
128+
mSeed = gRandom->Integer(conf.getStartSeed());
129+
} else {
130+
mSeed = seed;
131+
}
132+
LOG(info) << "Using seed: " << mSeed << " for HERWIG simulation";
133+
}
134+
135+
/// Destructor
136+
~GeneratorHerwig() = default;
137+
138+
/// Initialize the generator
139+
Bool_t Init() override
140+
{
141+
LOG(info) << "Initializing HERWIG7 Generator";
142+
if (mIsInputFile) {
143+
// Process .in file
144+
return initFromInputFile();
145+
} else {
146+
// Process .run file directly
147+
return initFromRunFile();
148+
}
149+
}
150+
151+
/// Generate event
152+
Bool_t generateEvent() override
153+
{
154+
if (!mEventGenerator) {
155+
LOG(error) << "Event generator not initialized";
156+
return kFALSE;
157+
}
158+
// Clear previous event particles
159+
hParticles.clear();
160+
161+
// Generate the event
162+
mPEGEvent = mEventGenerator->shoot();
163+
164+
if (!mPEGEvent)
165+
{
166+
LOG(error) << "Failed to generate event";
167+
return kFALSE;
168+
}
169+
170+
// Convert ThePEG event to TParticle format
171+
convertEvent(mPEGEvent);
172+
LOG(debug) << "Herwig7 generated " << hParticles.size() << " particles";
173+
174+
return kTRUE;
175+
}
176+
177+
/// Import particles for transport
178+
Bool_t importParticles() override
179+
{
180+
if (hParticles.empty()) {
181+
LOG(warning) << "No particles to import";
182+
return kFALSE;
183+
}
184+
185+
// Add particles to the primary generator
186+
for (const auto& particle : hParticles) {
187+
mParticles.push_back(particle);
188+
}
189+
190+
return kTRUE;
191+
}
192+
193+
private:
194+
std::string mConfigFile; ///< HERWIG config file (.in or .run)
195+
Bool_t mIsInputFile = false; ///< True for .in files, false for .run files
196+
ThePEG::EGPtr mEventGenerator; ///< ThePEG event generator
197+
std::vector<TParticle> hParticles; ///< Generated Herwig particles
198+
ThePEG::EventPtr mPEGEvent; ///< Pointer to Current event
199+
int mSeed = 0; ///< Random seed for Herwig
200+
201+
void printHerwigSearchPaths()
202+
{
203+
const auto &list = ThePEG::Repository::listReadDirs();
204+
205+
LOG(info) << "Append directories:\n";
206+
for (const auto &p : list)
207+
LOG(info) << " " << p << "\n";
208+
}
209+
210+
/// Initialize from .in file
211+
Bool_t initFromInputFile()
212+
{
213+
LOG(info) << "Initializing from .in file: " << mConfigFile;
214+
215+
using namespace ThePEG;
216+
SimpleHerwigUI ui(mConfigFile, Herwig::RunMode::READ, mSeed);
217+
Herwig::API::read(ui);
218+
// For debugging, print the search paths
219+
// printHerwigSearchPaths();
220+
// Currently the .run filename is set inside the .in file itself with
221+
// the line "saverun LHC EventGenerator" or similar. We assume this is the same as
222+
// the .in file name with .run extension, so change that string accordingly in your .in files.
223+
std::string runFile = mConfigFile;
224+
size_t pos = runFile.find_last_of('.');
225+
runFile.replace(pos, 4, ".run");
226+
pos = runFile.find_last_of('/');
227+
runFile = (pos != std::string::npos) ? runFile.substr(pos + 1) : runFile;
228+
mConfigFile = runFile;
229+
LOG(info) << "Generated run file: " << runFile;
230+
auto res = initFromRunFile();
231+
if (!res) {
232+
LOG(error) << "Failed to initialize from generated run file";
233+
return kFALSE;
234+
}
235+
return kTRUE;
236+
}
237+
238+
/// Initialize from .run file
239+
Bool_t initFromRunFile()
240+
{
241+
LOG(info) << "Initializing from .run file: " << mConfigFile;
242+
243+
using namespace ThePEG;
244+
245+
if (!std::ifstream(mConfigFile))
246+
{
247+
LOG(info) << "Run file does not exist: " << mConfigFile;
248+
return kFALSE;
249+
}
250+
SimpleHerwigUI runui(mConfigFile, Herwig::RunMode::RUN, mSeed);
251+
// Prepare the generator
252+
mEventGenerator = Herwig::API::prepareRun(runui);
253+
if (!mEventGenerator)
254+
{
255+
LOG(fatal) << "Error: prepareRun() returned null.";
256+
return kFALSE;
257+
}
258+
259+
mEventGenerator->initialize();
260+
LOG(info) << "Herwig generator initialized successfully.";
261+
return kTRUE;
262+
}
263+
264+
/// Convert ThePEG event to TParticle format
265+
void convertEvent(ThePEG::EventPtr event)
266+
{
267+
if (!event) return;
268+
269+
// Get all particles from the event
270+
const ThePEG::tPVector& particles = event->getFinalState();
271+
272+
for (size_t i = 0; i < particles.size(); ++i) {
273+
ThePEG::tPPtr particle = particles[i];
274+
if (!particle) continue;
275+
276+
// Get particle properties
277+
int pdgCode = particle->id();
278+
int status = getFinalStateStatus(particle);
279+
280+
// Get 4-momentum
281+
ThePEG::LorentzMomentum momentum = particle->momentum();
282+
double px = momentum.x() / ThePEG::GeV; // Convert to GeV
283+
double py = momentum.y() / ThePEG::GeV;
284+
double pz = momentum.z() / ThePEG::GeV;
285+
double e = momentum.e() / ThePEG::GeV;
286+
287+
// Get production vertex
288+
const ThePEG::LorentzPoint &vertex = particle->vertex();
289+
double vx = vertex.x() / ThePEG::mm; // Convert to mm
290+
double vy = vertex.y() / ThePEG::mm;
291+
double vz = vertex.z() / ThePEG::mm;
292+
double vt = vertex.t() / ThePEG::mm; // Convert to mm/c
293+
294+
// Create TParticle
295+
TParticle tparticle(
296+
pdgCode, status,
297+
-1, -1, -1, -1, // mother and daughter indices (to be set properly)
298+
px, py, pz, e,
299+
vx, vy, vz, vt
300+
);
301+
tparticle.SetStatusCode(o2::mcgenstatus::MCGenStatusEncoding(tparticle.GetStatusCode(), 0).fullEncoding);
302+
tparticle.SetBit(ParticleStatus::kToBeDone, //
303+
o2::mcgenstatus::getHepMCStatusCode(tparticle.GetStatusCode()) == 1);
304+
305+
hParticles.push_back(tparticle);
306+
}
307+
308+
LOG(debug) << "Converted " << hParticles.size() << " particles from HERWIG7 event";
309+
}
310+
311+
/// Determine final state status for particle
312+
int getFinalStateStatus(ThePEG::tPPtr particle)
313+
{
314+
// In HERWIG/ThePEG, check if particle is stable
315+
if (particle->children().empty()) {
316+
return 1; // Final state particle
317+
} else {
318+
return 2; // Intermediate particle
319+
}
320+
}
321+
};
322+
323+
} // namespace eventgen
324+
} // namespace o2
325+
326+
/// HERWIG7 generator from .in/.run file. If seed is -1, a random seed is chosen
327+
/// based on the SimConfig starting seed.
328+
FairGenerator* generateHerwig7(const std::string inputFile = "LHC.in", int seed = -1)
329+
{
330+
auto filePath = gSystem->ExpandPathName(inputFile.c_str());
331+
auto generator = new o2::eventgen::GeneratorHerwig(filePath, seed);
332+
return generator;
333+
}

‎MC/config/examples/ini/GeneratorHERWIG7.ini‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,3 @@
1+
[GeneratorExternal]
2+
fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/examples/herwig/generator_Herwig.C
3+
funcName=generateHerwig7("${O2DPG_MC_CONFIG_ROOT}/MC/config/examples/herwig/LHC.in")

0 commit comments

Comments
 (0)