Skip to content

Commit 8d1b770

Browse files
jackal1-66jackal1-66
committed
External Pythia8+Powheg generator
1 parent f503c23 commit 8d1b770

File tree

5 files changed

+288
-0
lines changed

5 files changed

+288
-0
lines changed

‎MC/config/PWGGAJE/external/generator/generator_pythia8_powheg.C‎

Lines changed: 102 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,102 @@
1+
R__ADD_INCLUDE_PATH($O2DPG_MC_CONFIG_ROOT)
2+
///#include "FairGenerator.h"
3+
//#include "Generators/GeneratorPythia8.h"
4+
#include "Pythia8/Pythia.h"
5+
#include "Generators/GeneratorPythia8Param.h"
6+
// Pythia8 generator with POWHEG
7+
//
8+
// Author: Marco Giacalone (marco.giacalone@cern.ch)
9+
10+
// o2-sim-dpl-eventgen --nEvents 10 --generator external\
11+
--configKeyValues "GeneratorExternal.fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/external/\
12+
generator/generator_pythia8_powheg.C;GeneratorExternal.funcName=\
13+
getGeneratorJEPythia8POWHEG(\"powheg.input\",\"${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/pythia8/generator/pythia8_powheg.cfg\")"
14+
// or with iniFile
15+
// o2-sim -g external --noGeant -n 2 -j 8 --configFile $O2DPG_MC_CONFIG_ROOT/MC/config/PWGGAJE/ini/GeneratorPythia8POWHEG.ini
16+
17+
namespace o2
18+
{
19+
namespace eventgen
20+
{
21+
22+
using namespace Pythia8;
23+
24+
class GeneratorJEPythia8POWHEG : public o2::eventgen::GeneratorPythia8
25+
{
26+
public:
27+
/// default constructor
28+
GeneratorJEPythia8POWHEG(std::string confpath = "pwgpath")
29+
{
30+
// Check if file exist and is not empty
31+
if (std::filesystem::exists(confpath) && std::filesystem::file_size(confpath) > 0) {
32+
// Copy the file to the current directory
33+
ifstream src(confpath);
34+
ofstream dst("powheg.input");
35+
gRandom->SetSeed(0);
36+
int seed = gRandom->Integer(900000000);
37+
bool isseed = false;
38+
bool isnumevts = false;
39+
std::string line;
40+
while (std::getline(src, line)) {
41+
if (line.find("iseed") != std::string::npos)
42+
{
43+
// Set the seed to the random number
44+
line = "iseed " + std::to_string(seed);
45+
isseed = true;
46+
}
47+
if (line.find("numevts") != std::string::npos)
48+
{
49+
// Set the number of events to the number of events defined in the configuration
50+
line = "numevts " + std::to_string(mSimConfig.getNEvents());
51+
// replace it in the file
52+
isnumevts = true;
53+
}
54+
dst << line << std::endl;
55+
}
56+
if (!isseed) {
57+
dst << "iseed " << seed << std::endl;
58+
}
59+
if (!isnumevts) {
60+
dst << "numevts " << mSimConfig.getNEvents() << std::endl;
61+
}
62+
src.close();
63+
dst.close();
64+
} else {
65+
LOG(fatal) << "POWHEG configuration file not found or empty" << std::endl;
66+
exit(1);
67+
}
68+
// Generate the POWHEG events
69+
std::string cmd = "pwhg_main_hvq";
70+
system(cmd.c_str());
71+
};
72+
73+
private:
74+
o2::conf::SimConfig mSimConfig = o2::conf::SimConfig::Instance(); // local sim config object
75+
};
76+
77+
} // namespace eventgen
78+
} // namespace o2
79+
80+
/** generator instance and settings **/
81+
82+
FairGenerator *getGeneratorJEPythia8POWHEG(std::string powhegconf = "pwgpath", std::string pythia8conf = "")
83+
{
84+
using namespace o2::eventgen;
85+
// Check if the conf path contains the O2DPG configuration environment variable
86+
if (powhegconf.find("${O2DPG_MC_CONFIG_ROOT}") != std::string::npos) {
87+
powhegconf.replace(powhegconf.find("${O2DPG_MC_CONFIG_ROOT}"), 23, std::getenv("O2DPG_MC_CONFIG_ROOT"));
88+
}
89+
auto myGen = new GeneratorJEPythia8POWHEG(powhegconf);
90+
if(GeneratorPythia8Param::Instance().config.empty() && pythia8conf.empty()) {
91+
LOG(fatal) << "No configuration provided for Pythia8";
92+
}
93+
else if (!pythia8conf.empty())
94+
{
95+
// Force the configuration for Pythia8 in case it is provided.
96+
// Useful for setting up the generator in the hybrid configuration
97+
// making it more versatile and not relying entirely on the parameters provided
98+
// by ini file or static parameters
99+
myGen->setConfig(pythia8conf);
100+
}
101+
return myGen;
102+
}

‎MC/config/PWGGAJE/external/powheg.input‎

Lines changed: 52 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,52 @@
1+
!randomseed 352345 ! uncomment to set the random seed to a value of your choice.
2+
! It generates the call RM48IN(352345,0,0) (see the RM48 manual).
3+
! THIS MAY ONLY AFFECTS THE GENERATION OF POWHEG EVENTS!
4+
! If POWHEG is interfaced to a shower MC, refer to the shower MC
5+
! documentation to set its seed.
6+
7+
!Heavy flavour production parameters
8+
9+
ih1 1 ! hadron 1
10+
ih2 1 ! hadron 2
11+
#ndns1 131 ! pdf for hadron 1 (hvqpdf numbering)
12+
#ndns2 131 ! pdf for hadron 2
13+
!lhans1 11000 ! pdf set for hadron 1 (LHA numbering)
14+
!lhans2 11000 ! pdf set for hadron 2 (LHA numbering)
15+
!ebeam1 3500 ! energy of beam 1
16+
!ebeam2 3500 ! energy of beam 2
17+
!qmass 1.5 ! mass of heavy quark in GeV
18+
!facscfact 1 ! factorization scale factor: mufact=muref*facscfact
19+
!renscfact 1 ! renormalization scale factor: muren=muref*renscfact
20+
#fixedscale 1 ! use ref. scale=qmass (default 0, use running scale)
21+
22+
! Parameters to allow-disallow use of stored data
23+
use-old-grid 1 ! if 1 use old grid if file pwggrids.dat is present (# 1: regenerate)
24+
use-old-ubound 1 ! if 1 use norm of upper bounding function stored in pwgubound.dat, if present; # 1: regenerate
25+
26+
!ncall1 50000 ! number of calls for initializing the integration grid
27+
!itmx1 5 ! number of iterations for initializing the integration grid
28+
!ncall2 100000 ! number of calls for computing the integral and finding upper bound
29+
!itmx2 5 ! number of iterations for computing the integral and finding upper bound
30+
foldcsi 5 ! number of folds on x integration
31+
foldy 5 ! number of folds on y integration
32+
foldphi 1 ! number of folds on phi integration
33+
nubound 500000 ! number of bbarra calls to setup norm of upper bounding function
34+
iymax 1 ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
35+
!ixmax 1 ! <= 10, normalization of upper bounding function in iunorm X iunorm square in y, log(m2qq)
36+
xupbound 2 ! increase upper bound for radiation generation
37+
qmass 1.5
38+
ncall1 50000
39+
itmx1 5
40+
ncall2 100000
41+
itmx2 5
42+
facscfact 1
43+
renscfact 1
44+
bornktmin 0
45+
bornsuppfact 0
46+
storemintupb 0
47+
lhans1 11000
48+
lhans2 11000
49+
ebeam1 6800
50+
ebeam2 6800
51+
iseed 389759591
52+
numevts 4

‎MC/config/PWGGAJE/ini/GeneratorPythia8POWHEG.ini‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
### The external generator derives from GeneratorPythia8.
2+
## The generator allows to run Pythia8 with POWHEG
3+
[GeneratorExternal]
4+
fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/external/generator/generator_pythia8_powheg.C
5+
funcName=getGeneratorJEPythia8POWHEG("${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/external/powheg.input")
6+
7+
[GeneratorPythia8]
8+
config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGGAJE/pythia8/generator/pythia8_powheg.cfg
9+

‎MC/config/PWGGAJE/ini/tests/GeneratorPythia8POWHEG.C‎

Lines changed: 46 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,46 @@
1+
int External() {
2+
std::string path{"o2sim_Kine.root"};
3+
std::ifstream powhegconf("powheg.input");
4+
if (!powhegconf) {
5+
std::cerr << "POWHEG configuration file not found\n";
6+
return 1;
7+
}
8+
std::ifstream powhegout("pwgevents.lhe");
9+
if (!powhegout) {
10+
std::cerr << "POWHEG output file not found\n";
11+
return 1;
12+
}
13+
powhegout.close();
14+
15+
TFile file(path.c_str(), "READ");
16+
if (file.IsZombie()) {
17+
std::cerr << "Cannot open ROOT file " << path << "\n";
18+
return 1;
19+
}
20+
21+
auto tree = (TTree *)file.Get("o2sim");
22+
auto nEvents = tree->GetEntries();
23+
std::string line;
24+
int nevpowheg = -1;
25+
while (std::getline(powhegconf, line)) {
26+
if (line.find("numevts") != std::string::npos) {
27+
// Read the number right after numevts
28+
auto pos = line.find("numevts");
29+
nevpowheg = std::stoi(line.substr(pos + 7));
30+
if (nevpowheg != nEvents) {
31+
std::cerr << "Number of events in POWHEG configuration file " << nevpowheg
32+
<< " does not match the simulated number of events "
33+
<< nEvents << "\n";
34+
return 1;
35+
}
36+
}
37+
}
38+
if (nevpowheg == -1) {
39+
std::cerr << "Number of events not found in POWHEG configuration file\n";
40+
return 1;
41+
}
42+
powhegconf.close();
43+
file.Close();
44+
45+
return 0;
46+
}

‎MC/config/PWGGAJE/pythia8/generator/pythia8_powheg.cfg‎

Lines changed: 79 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,79 @@
1+
Beams:eCM 13000. # GeV
2+
### processes
3+
Next:numberShowLHA = 1
4+
Next:numberShowInfo = 1
5+
Next:numberShowProcess = 1
6+
Next:numberShowEvent = 1
7+
Main:timesAllowErrors = 10
8+
! Read LHE file from POWHEG
9+
Beams:frameType = 4
10+
Beams:LHEF = pwgevents.lhe
11+
12+
! Number of outgoing particles of POWHEG Born level process
13+
! (i.e. not counting additional POWHEG radiation)
14+
POWHEG:nFinal = 2
15+
16+
! How vetoing is performed:
17+
! 0 - No vetoing is performed (userhooks are not loaded)
18+
! 1 - Showers are started at the kinematical limit.
19+
! Emissions are vetoed if pTemt > pThard.
20+
! See also POWHEG:vetoCount below
21+
POWHEG:veto = 1
22+
23+
! After 'vetoCount' accepted emissions in a row, no more emissions
24+
! are checked. 'vetoCount = 0' means that no emissions are checked.
25+
! Use a very large value, e.g. 10000, to have all emissions checked.
26+
POWHEG:vetoCount = 3
27+
28+
! Selection of pThard (note, for events where there is no
29+
! radiation, pThard is always set to be SCALUP):
30+
! 0 - pThard = SCALUP (of the LHA/LHEF standard)
31+
! 1 - the pT of the POWHEG emission is tested against all other
32+
! incoming and outgoing partons, with the minimal value chosen
33+
! 2 - the pT of all final-state partons is tested against all other
34+
! incoming and outgoing partons, with the minimal value chosen
35+
POWHEG:pThard = 2
36+
37+
! Selection of pTemt:
38+
! 0 - pTemt is pT of the emitted parton w.r.t. radiating parton
39+
! 1 - pT of the emission is checked against all incoming and outgoing
40+
! partons. pTemt is set to the minimum of these values
41+
! 2 - the pT of all final-state partons is tested against all other
42+
! incoming and outgoing partons, with the minimal value chosen
43+
! WARNING: the choice here can give significant variations in the final
44+
! distributions, notably in the tail to large pT values.
45+
POWHEG:pTemt = 0
46+
47+
! Selection of emitted parton for FSR
48+
! 0 - Pythia definition of emitted
49+
! 1 - Pythia definition of radiator
50+
! 2 - Random selection of emitted or radiator
51+
! 3 - Both are emitted and radiator are tried
52+
POWHEG:emitted = 0
53+
54+
! pT definitions
55+
! 0 - POWHEG ISR pT definition is used for both ISR and FSR
56+
! 1 - POWHEG ISR pT and FSR d_ij definitions
57+
! 2 - Pythia definitions
58+
POWHEG:pTdef = 1
59+
60+
! MPI vetoing
61+
! 0 - No MPI vetoing is done
62+
! 1 - When there is no radiation, MPIs with a scale above pT_1 are vetoed,
63+
! else MPIs with a scale above (pT_1 + pT_2 + pT_3) / 2 are vetoed
64+
POWHEG:MPIveto = 0
65+
! Note that POWHEG:MPIveto = 1 should be combined with
66+
! MultipartonInteractions:pTmaxMatch = 2
67+
! which here is taken care of in main31.cc.
68+
69+
! QED vetoing
70+
! 0 - No QED vetoing is done for pTemt > 0.
71+
! 1 - QED vetoing is done for pTemt > 0.
72+
! 2 - QED vetoing is done for pTemt > 0. If a photon is found
73+
! with pT>pThard from the Born level process, the event is accepted
74+
! and no further veto of this event is allowed (for any pTemt).
75+
POWHEG:QEDveto = 2
76+
77+
### decays
78+
ParticleDecays:limitTau0 on
79+
ParticleDecays:tau0Max 10.

0 commit comments

Comments
 (0)