Page MenuHomeHEPForge
Files F8725696

No OneTemporary
Actions

View Options

diff --git a/FixedOrderGen/include/Beam.hh b/FixedOrderGen/include/Beam.hh
new file mode 100644
index 0000000..48a0f47
--- /dev/null
+++ b/FixedOrderGen/include/Beam.hh
@@ -0,0 +1,14 @@
+#pragma once
+
+#include <array>
+
+#include "RHEJ/PDG_codes.hh"
+
+namespace HEJFOG{
+ struct Beam{
+ double energy;
+ std::array<RHEJ::ParticleID, 2> particles = {
+ RHEJ::pid::proton, RHEJ::pid::proton
+ };
+ };
+}
diff --git a/FixedOrderGen/include/EventGenerator.hh b/FixedOrderGen/include/EventGenerator.hh
new file mode 100644
index 0000000..29edd01
--- /dev/null
+++ b/FixedOrderGen/include/EventGenerator.hh
@@ -0,0 +1,46 @@
+#pragma once
+
+#include "RHEJ/PDF.hh"
+#include "RHEJ/MatrixElement.hh"
+
+#include "JetParameters.hh"
+#include "Process.hh"
+#include "Beam.hh"
+
+namespace RHEJ{
+ class Event;
+ class HiggsCouplingSettings;
+ class ScaleGenerator;
+}
+
+namespace HEJFOG{
+ class EventGenerator{
+ public:
+ EventGenerator(
+ Process process,
+ Beam beam,
+ RHEJ::ScaleGenerator & scale_gen,
+ JetParameters jets,
+ int pdf_id,
+ bool unordered,
+ RHEJ::HiggsCouplingSettings Higgs_coupling
+ );
+
+ RHEJ::Event gen_event();
+
+ int status() const {
+ return status_;
+ }
+
+ private:
+ RHEJ::PDF pdf_;
+ RHEJ::MatrixElement ME_;
+ RHEJ::ScaleGenerator & scale_gen_;
+ Process process_;
+ JetParameters jets_;
+ Beam beam_;
+ bool unordered_;
+ int status_;
+ };
+
+}
diff --git a/FixedOrderGen/include/JetParameters.hh b/FixedOrderGen/include/JetParameters.hh
new file mode 100644
index 0000000..81e7354
--- /dev/null
+++ b/FixedOrderGen/include/JetParameters.hh
@@ -0,0 +1,9 @@
+#pragma once
+
+namespace HEJFOG{
+ struct JetParameters{
+ fastjet::JetDefinition def;
+ double min_pt;
+ double max_y;
+ };
+}
diff --git a/FixedOrderGen/include/config.hh b/FixedOrderGen/include/config.hh
index c52947a..1c9f0f4 100644
--- a/FixedOrderGen/include/config.hh
+++ b/FixedOrderGen/include/config.hh
@@ -1,47 +1,34 @@
#pragma once
-#include <array>
-
#include "yaml-cpp/yaml.h"
#include "fastjet/JetDefinition.hh"
#include "RHEJ/HiggsCouplingSettings.hh"
-#include "RHEJ/PDG_codes.hh"
#include "RHEJ/optional.hh"
#include "RHEJ/config.hh"
#include "RHEJ/output_formats.hh"
#include "RHEJ/exceptions.hh"
#include "Process.hh"
+#include "JetParameters.hh"
+#include "Beam.hh"
namespace HEJFOG{
- struct JetParameters{
- fastjet::JetDefinition def;
- double min_pt;
- double max_y;
- };
-
- struct Beam{
- double energy;
- std::array<RHEJ::ParticleID, 2> particles = {
- RHEJ::pid::proton, RHEJ::pid::proton
- };
- };
struct Config{
Process process;
int trials;
JetParameters jets;
Beam beam;
int pdf_id;
bool unordered = true;
YAML::Node analysis_parameters;
RHEJ::ScaleGenerator scale_gen;
std::vector<RHEJ::OutputFile> output;
RHEJ::optional<std::string> RanLux_init;
RHEJ::HiggsCouplingSettings Higgs_coupling;
};
Config load_config(std::string const & config_file);
}
diff --git a/FixedOrderGen/src/EventGenerator.cc b/FixedOrderGen/src/EventGenerator.cc
new file mode 100644
index 0000000..187854b
--- /dev/null
+++ b/FixedOrderGen/src/EventGenerator.cc
@@ -0,0 +1,72 @@
+#include "EventGenerator.hh"
+
+#include "Process.hh"
+#include "Beam.hh"
+#include "JetParameters.hh"
+#include "PhaseSpacePoint.hh"
+
+#include "RHEJ/Event.hh"
+#include "RHEJ/config.hh"
+
+namespace HEJFOG{
+ EventGenerator::EventGenerator(
+ Process process,
+ Beam beam,
+ RHEJ::ScaleGenerator & scale_gen,
+ JetParameters jets,
+ int pdf_id,
+ bool unordered,
+ RHEJ::HiggsCouplingSettings Higgs_coupling
+ ):
+ pdf_{pdf_id, beam.particles[0], beam.particles[1]},
+ ME_{
+ [this](double mu){ return pdf_.Halphas(mu); },
+ jets.def, jets.min_pt,
+ false,
+ std::move(Higgs_coupling)
+ },
+ scale_gen_{scale_gen},
+ process_{std::move(process)},
+ jets_{std::move(jets)},
+ beam_{std::move(beam)},
+ unordered_{unordered}
+ {
+ }
+
+ RHEJ::Event EventGenerator::gen_event(){
+ HEJFOG::PhaseSpacePoint psp{
+ process_,
+ jets_.def, jets_.min_pt, jets_.max_y,
+ pdf_, beam_.energy
+ };
+ status_ = psp.status();
+ if(status_ != 0) return {};
+
+ RHEJ::Event ev = scale_gen_(
+ RHEJ::Event{
+ to_UnclusteredEvent(std::move(psp)),
+ jets_.def, jets_.min_pt
+ }
+ );
+
+ const double shat = RHEJ::shat(ev);
+ const double xa = (ev.incoming()[0].E()-ev.incoming()[0].pz())/(2.*beam_.energy);
+ const double xb = (ev.incoming()[1].E()+ev.incoming()[1].pz())/(2.*beam_.energy);
+
+ // evaluate matrix element on this point
+ ev.central().weight *= ME_.tree(
+ ev.central().mur, ev.incoming(), ev.outgoing(), false
+ )/(shat*shat);
+ ev.central().weight *= pdf_.pdfpt(0,xa,ev.central().muf, ev.incoming()[0].type);
+ ev.central().weight *= pdf_.pdfpt(0,xb,ev.central().muf, ev.incoming()[1].type);
+ for(auto & var: ev.variations()){
+ var.weight *= ME_.tree(
+ var.mur, ev.incoming(), ev.outgoing(), false
+ )/(shat*shat);
+ var.weight *= pdf_.pdfpt(0,xa,var.muf, ev.incoming()[0].type);
+ var.weight *= pdf_.pdfpt(0,xb,var.muf, ev.incoming()[1].type);
+ }
+ return ev;
+ }
+
+}
diff --git a/FixedOrderGen/src/main.cc b/FixedOrderGen/src/main.cc
index c954889..47b1c64 100644
--- a/FixedOrderGen/src/main.cc
+++ b/FixedOrderGen/src/main.cc
@@ -1,195 +1,165 @@
/**
* Name: main.cc
* Authors: Jeppe R. Andersen
*/
#include <fstream>
#include <algorithm>
#include <memory>
#include <chrono>
#include <iostream>
#include <TROOT.h>
#include <TFile.h>
#include <TH1.h>
#include "yaml-cpp/yaml.h"
#include "config.hh"
#include "LHEF/LHEF.h"
#include "RHEJ/CombinedEventWriter.hh"
#include "RHEJ/get_analysis.hh"
#include "RHEJ/utility.hh"
-#include "RHEJ/PDF.hh"
//#include "RHEJ/EventReweighter.hh"
#include "RHEJ/stream.hh"
-#include "RHEJ/MatrixElement.hh"
#include "RHEJ/LesHouchesWriter.hh"
+
+#include "EventGenerator.hh"
#include "PhaseSpacePoint.hh"
namespace{
constexpr auto banner =
" __ ___ __ ______ __ __ \n"
" / / / (_)___ _/ /_ / ____/___ ___ _________ ___ __ / /__ / /______ \n"
" / /_/ / / __ `/ __ \\ / __/ / __ \\/ _ \\/ ___/ __ `/ / / / __ / / _ \\/ __/ ___/ \n"
" / __ / / /_/ / / / / / /___/ / / / __/ / / /_/ / /_/ / / /_/ / __/ /_(__ ) \n"
" /_/ /_/_/\\__, /_/ /_/ /_____/_/ /_/\\___/_/ \\__, /\\__, / \\____/\\___/\\__/____/ \n"
" ____///__/ __ ____ ///__//____/ ______ __ \n"
" / ____(_) _____ ____/ / / __ \\_________/ /__ _____ / ____/__ ____ ___ _________ _/ /_____ _____\n"
" / /_ / / |/_/ _ \\/ __ / / / / / ___/ __ / _ \\/ ___/ / / __/ _ \\/ __ \\/ _ \\/ ___/ __ `/ __/ __ \\/ ___/\n"
" / __/ / /> </ __/ /_/ / / /_/ / / / /_/ / __/ / / /_/ / __/ / / / __/ / / /_/ / /_/ /_/ / / \n"
" /_/ /_/_/|_|\\___/\\__,_/ \\____/_/ \\__,_/\\___/_/ \\____/\\___/_/ /_/\\___/_/ \\__,_/\\__/\\____/_/ \n"
;
constexpr double invGeV2_to_pb = 389379292.;
}
HEJFOG::Config load_config(char const * filename){
try{
return HEJFOG::load_config(filename);
}
catch(std::exception const & exc){
std::cerr << "Error: " << exc.what() << '\n';
std::exit(EXIT_FAILURE);
}
}
std::unique_ptr<RHEJ::Analysis> get_analysis(
YAML::Node const & parameters
){
try{
return RHEJ::get_analysis(parameters);
}
catch(std::exception const & exc){
std::cerr << "Failed to load analysis: " << exc.what() << '\n';
std::exit(EXIT_FAILURE);
}
}
std::string progress_bar(){
std::string result = "0% ";
for(int i = 10; i <= 100; i+= 10){
result += " " + std::to_string(i) + "%";
}
result += "\n|";
for(int i = 10; i <= 100; i+= 10){
result += "---------|";
}
return result + '\n';
}
int main(int argn, char** argv) {
if (argn < 2) {
std::cerr << "\n# Usage:\n.FOgen config_file\n";
return EXIT_FAILURE;
}
std::cout << banner;
fastjet::ClusterSequence::print_banner();
using clock = std::chrono::system_clock;
const auto start_time = clock::now();
// read configuration
- const auto config = load_config(argv[1]);
+ auto config = load_config(argv[1]);
std::unique_ptr<RHEJ::Analysis> analysis = get_analysis(
config.analysis_parameters
);
assert(analysis != nullptr);
- RHEJ::PDF pdf{
+ HEJFOG::EventGenerator generator{
+ config.process,
+ config.beam,
+ config.scale_gen,
+ config.jets,
config.pdf_id,
- config.beam.particles[0], config.beam.particles[1]
- };
-
- // Generate a matrix element:
- RHEJ::MatrixElement ME{
- [&pdf](double mu){ return pdf.Halphas(mu); },
- config.jets.def, config.jets.min_pt,
- false,
+ config.unordered,
config.Higgs_coupling
};
//TODO: fix Les Houches init block (HEPRUP)
LHEF::HEPRUP lhefheprup;
lhefheprup.IDBMUP=std::pair<long,long>(config.beam.particles[0], config.beam.particles[1]);
lhefheprup.EBMUP=std::make_pair(config.beam.energy, config.beam.energy);
lhefheprup.PDFGUP=std::make_pair(0,0);
lhefheprup.PDFSUP=std::make_pair(config.pdf_id,config.pdf_id);
lhefheprup.NPRUP=1;
lhefheprup.XSECUP=std::vector<double>(1.);
lhefheprup.XERRUP=std::vector<double>(1.);
lhefheprup.LPRUP=std::vector<int>{1};
RHEJ::CombinedEventWriter writer{config.output, lhefheprup};
if(config.RanLux_init){
HEJFOG::PhaseSpacePoint::reset_ranlux(*config.RanLux_init);
}
// Book root histogram for status
TFile hfile("GenStatus.root","RECREATE","Generator status");
TROOT simple("FOG","Output from HEJFOG");
TH1D *hmstatus;
hmstatus = new TH1D("mstatus","mstatus",25,-.5,24.5);
// Perform N trial generations
int iprint = 0;
std::cout << '\n' << progress_bar();
for (int trials = 0; trials< config.trials;trials++){
if (trials==iprint) {
std::cout << ".";
std::cout.flush();
iprint+=(int) config.trials/100;
}
- // Generate phase space point
- HEJFOG::PhaseSpacePoint psp{
- config.process,
- config.jets.def,config.jets.min_pt, config.jets.max_y,
- pdf, config.beam.energy
- };
- hmstatus->Fill(psp.status());
- if (psp.status()!=0) continue;
-
- RHEJ::Event ev = config.scale_gen(
- RHEJ::Event{
- to_UnclusteredEvent(std::move(psp)),
- config.jets.def, config.jets.min_pt
- }
- );
- const double shat = RHEJ::shat(ev);
- const double xa = (ev.incoming()[0].E()-ev.incoming()[0].pz())/(2.*config.beam.energy);
- const double xb = (ev.incoming()[1].E()+ev.incoming()[1].pz())/(2.*config.beam.energy);
-
- // evaluate matrix element on this point
- ev.central().weight *= ME.tree(
- ev.central().mur, ev.incoming(), ev.outgoing(), false
- )/(shat*shat);
- ev.central().weight *= pdf.pdfpt(0,xa,ev.central().muf, ev.incoming()[0].type);
- ev.central().weight *= pdf.pdfpt(0,xb,ev.central().muf, ev.incoming()[1].type);
- ev.central().weight *= invGeV2_to_pb;
+ auto ev = generator.gen_event();
+ hmstatus->Fill(generator.status());
+ if(generator.status() != 0) continue;
ev.central().weight /= config.trials;
+ ev.central().weight *= invGeV2_to_pb;
for(auto & var: ev.variations()){
- var.weight *= ME.tree(
- var.mur, ev.incoming(), ev.outgoing(), false
- )/(shat*shat);
- var.weight *= pdf.pdfpt(0,xa,var.muf, ev.incoming()[0].type);
- var.weight *= pdf.pdfpt(0,xb,var.muf, ev.incoming()[1].type);
- var.weight *= invGeV2_to_pb;
var.weight /= config.trials;
+ var.weight *= invGeV2_to_pb;
}
+
if(analysis->pass_cuts(ev)){
analysis->fill(ev);
writer.write(ev);
}
-
} // main event loop
std::cout << std::endl;
std::chrono::duration<double> run_time = (clock::now() - start_time);
std::cout << "\nTask Runtime: " << run_time.count() << " seconds.\n";
hfile.Write();
hfile.Close();
}
diff --git a/FixedOrderGen/t/h_2j.cc b/FixedOrderGen/t/h_2j.cc
index c2993fb..847f792 100644
--- a/FixedOrderGen/t/h_2j.cc
+++ b/FixedOrderGen/t/h_2j.cc
@@ -1,79 +1,56 @@
#ifdef NDEBUG
#undef NDEBUG
#endif
#include <algorithm>
#include <cmath>
#include <cassert>
#include <iostream>
#include "config.hh"
-#include "PhaseSpacePoint.hh"
+#include "EventGenerator.hh"
#include "RHEJ/PDF.hh"
#include "RHEJ/MatrixElement.hh"
using namespace HEJFOG;
int main(){
constexpr double invGeV2_to_pb = 389379292.;
constexpr double xs_ref = 2.0304; //calculated with "old" HEJ svn r3364
- const auto config = load_config("config_h_2j.yml");
+ auto config = load_config("config_h_2j.yml");
- RHEJ::PDF pdf{
+ HEJFOG::EventGenerator generator{
+ config.process,
+ config.beam,
+ config.scale_gen,
+ config.jets,
config.pdf_id,
- config.beam.particles[0], config.beam.particles[1]
- };
-
- RHEJ::MatrixElement ME{
- [&pdf](double mu){ return pdf.Halphas(mu); },
- config.jets.def, config.jets.min_pt,
- false,
+ config.unordered,
config.Higgs_coupling
};
double xs = 0., xs_err = 0.;
for (int trials = 0; trials < config.trials; ++trials){
- HEJFOG::PhaseSpacePoint psp{
- config.process,
- config.jets.def,config.jets.min_pt, config.jets.max_y,
- pdf, config.beam.energy
- };
- if (psp.status()!=0) continue;
-
- RHEJ::Event ev = config.scale_gen(
- RHEJ::Event{
- to_UnclusteredEvent(std::move(psp)),
- config.jets.def, config.jets.min_pt
- }
- );
- const double shat = RHEJ::shat(ev);
- const double xa = (ev.incoming()[0].E()-ev.incoming()[0].pz())/(2.*config.beam.energy);
- const double xb = (ev.incoming()[1].E()+ev.incoming()[1].pz())/(2.*config.beam.energy);
-
- // evaluate matrix element on this point
- ev.central().weight *= ME.tree(
- ev.central().mur, ev.incoming(), ev.outgoing(), false
- )/(shat*shat);
- ev.central().weight *= pdf.pdfpt(0,xa,ev.central().muf, ev.incoming()[0].type);
- ev.central().weight *= pdf.pdfpt(0,xb,ev.central().muf, ev.incoming()[1].type);
+ auto ev = generator.gen_event();
+ if(generator.status() != 0) continue;
ev.central().weight *= invGeV2_to_pb;
ev.central().weight /= config.trials;
const auto the_Higgs = std::find_if(
begin(ev.outgoing()), end(ev.outgoing()),
[](RHEJ::Sparticle const & p){ return p.type == RHEJ::ParticleID::h; }
);
assert(the_Higgs != end(ev.outgoing()));
if(std::abs(the_Higgs->rapidity()) > 5.) continue;
xs += ev.central().weight;
xs_err += ev.central().weight*ev.central().weight;
}
xs_err = std::sqrt(xs_err);
std::cout << xs_ref << " ~ " << xs << " +- " << xs_err << '\n';
assert(std::abs(xs - xs_ref) < 3*xs_err);
assert(xs_err < 0.01*xs);
}
diff --git a/FixedOrderGen/t/h_3j.cc b/FixedOrderGen/t/h_3j.cc
index e169fbf..e0b7c48 100644
--- a/FixedOrderGen/t/h_3j.cc
+++ b/FixedOrderGen/t/h_3j.cc
@@ -1,80 +1,57 @@
#ifdef NDEBUG
#undef NDEBUG
#endif
#include <algorithm>
#include <cmath>
#include <cassert>
#include <iostream>
#include "config.hh"
-#include "PhaseSpacePoint.hh"
+#include "EventGenerator.hh"
#include "RHEJ/PDF.hh"
#include "RHEJ/MatrixElement.hh"
using namespace HEJFOG;
int main(){
constexpr double invGeV2_to_pb = 389379292.;
constexpr double xs_ref = 1.0888; //calculated with "old" HEJ svn r3364
auto config = load_config("config_h_2j.yml");
config.process.njets = 3;
- RHEJ::PDF pdf{
+ HEJFOG::EventGenerator generator{
+ config.process,
+ config.beam,
+ config.scale_gen,
+ config.jets,
config.pdf_id,
- config.beam.particles[0], config.beam.particles[1]
- };
-
- RHEJ::MatrixElement ME{
- [&pdf](double mu){ return pdf.Halphas(mu); },
- config.jets.def, config.jets.min_pt,
- false,
+ config.unordered,
config.Higgs_coupling
};
double xs = 0., xs_err = 0.;
for (int trials = 0; trials < config.trials; ++trials){
- HEJFOG::PhaseSpacePoint psp{
- config.process,
- config.jets.def,config.jets.min_pt, config.jets.max_y,
- pdf, config.beam.energy
- };
- if (psp.status()!=0) continue;
-
- RHEJ::Event ev = config.scale_gen(
- RHEJ::Event{
- to_UnclusteredEvent(std::move(psp)),
- config.jets.def, config.jets.min_pt
- }
- );
- const double shat = RHEJ::shat(ev);
- const double xa = (ev.incoming()[0].E()-ev.incoming()[0].pz())/(2.*config.beam.energy);
- const double xb = (ev.incoming()[1].E()+ev.incoming()[1].pz())/(2.*config.beam.energy);
-
- // evaluate matrix element on this point
- ev.central().weight *= ME.tree(
- ev.central().mur, ev.incoming(), ev.outgoing(), false
- )/(shat*shat);
- ev.central().weight *= pdf.pdfpt(0,xa,ev.central().muf, ev.incoming()[0].type);
- ev.central().weight *= pdf.pdfpt(0,xb,ev.central().muf, ev.incoming()[1].type);
+ auto ev = generator.gen_event();
+ if(generator.status() != 0) continue;
ev.central().weight *= invGeV2_to_pb;
ev.central().weight /= config.trials;
const auto the_Higgs = std::find_if(
begin(ev.outgoing()), end(ev.outgoing()),
[](RHEJ::Sparticle const & p){ return p.type == RHEJ::ParticleID::h; }
);
assert(the_Higgs != end(ev.outgoing()));
if(std::abs(the_Higgs->rapidity()) > 5.) continue;
xs += ev.central().weight;
xs_err += ev.central().weight*ev.central().weight;
}
xs_err = std::sqrt(xs_err);
std::cout << xs_ref << " ~ " << xs << " +- " << xs_err << '\n';
assert(std::abs(xs - xs_ref) < 3*xs_err);
assert(xs_err < 0.02*xs);
}

File Metadata

Mime Type
text/x-diff
Expires
Tue, Jan 21, 2:18 AM (1 d, 18 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4243648
Default Alt Text
(17 KB)

Event Timeline