diff --git a/src/EventReader.cc b/src/EventReader.cc
index 08d65f7..5015b74 100644
--- a/src/EventReader.cc
+++ b/src/EventReader.cc
@@ -1,64 +1,34 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/EventReader.hh"
#include <iostream>
#include <stdexcept>
#include "LHEF/LHEF.h"
#include "HEJ/ConfigFlags.hh"
#include "HEJ/HDF5Reader.hh"
#include "HEJ/LesHouchesReader.hh"
#include "HEJ/Version.hh"
-// namespace {
-
-// HEJ::generator get_generator(
-// LHEF::HEPRUP const & heprup, std::string const & header
-// ){
-// std::cout<<"header: "<<header<<std::endl;
-// // try getting generator name from specific tag
-// if(!heprup.generators.empty()){
-// std::string const & gen_name = heprup.generators.back().name;
-// if(gen_name == "HEJ" || gen_name == HEJ::Version::String())
-// return HEJ::generator::HEJ;
-// if(gen_name == "HEJ Fixed Order Generation")
-// return HEJ::generator::HEJFOG;
-// if(gen_name == "MadGraph5_aMC@NLO") return HEJ::generator::MG;
-// std::cerr << "Unknown LHE Generator " << gen_name
-// << " using default LHE interface.\n";
-// return HEJ::generator::unknown;
-// }
-// // The generator tag is not always used -> check by hand
-// if(header.find("generated with HEJ")!=std::string::npos)
-// return HEJ::generator::HEJ;
-// if(header.find("# created by SHERPA")!=std::string::npos)
-// return HEJ::generator::Sherpa;
-// if(header.find("<MGVersion>")!=std::string::npos)
-// return HEJ::generator::MG;
-// std::cerr<<"Could not determine LHE Generator using default LHE interface.\n";
-// return HEJ::generator::unknown;
-// }
-// } // namespace
-
namespace HEJ {
std::unique_ptr<EventReader> make_reader(std::string const & filename) {
try {
auto reader{ std::make_unique<HEJLHEReader>(filename) };
return reader;
}
catch(std::runtime_error&) {
#ifdef HEJ_BUILD_WITH_HDF5
return std::make_unique<HDF5Reader>(filename);
#else
throw;
#endif
}
}
} // namespace HEJ
diff --git a/src/LesHouchesReader.cc b/src/LesHouchesReader.cc
index 5f730af..57b767e 100644
--- a/src/LesHouchesReader.cc
+++ b/src/LesHouchesReader.cc
@@ -1,112 +1,82 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2022
* \copyright GPLv2 or later
*/
#include "HEJ/LesHouchesReader.hh"
#include "HEJ/Version.hh"
#include <vector>
namespace HEJ {
HEJ::generator get_generator(
LHEF::HEPRUP const & heprup, std::string const & header
){
// try getting generator name from specific tag
if(!heprup.generators.empty()){
std::string const & gen_name = heprup.generators.back().name;
if(gen_name == "HEJ" || gen_name == HEJ::Version::String())
return HEJ::generator::HEJ;
if(gen_name == "HEJ Fixed Order Generation")
return HEJ::generator::HEJFOG;
if(gen_name == "MadGraph5_aMC@NLO") return HEJ::generator::MG;
std::cerr << "Unknown LHE Generator " << gen_name
<< " using default LHE interface.\n";
return HEJ::generator::unknown;
}
// The generator tag is not always used -> check by hand
if(header.find("generated with HEJ")!=std::string::npos)
return HEJ::generator::HEJ;
if(header.find("# created by SHERPA")!=std::string::npos)
return HEJ::generator::Sherpa;
if(header.find("<MGVersion>")!=std::string::npos)
return HEJ::generator::MG;
std::cerr<<"Could not determine LHE Generator using default LHE interface.\n";
return HEJ::generator::unknown;
}
HEJLHEReader::HEJLHEReader(std::string const & filename):
LesHouchesReader{filename},
num_trials_(10,0),
num_events_(0),
CalculateXSECUP{true}
{
- // LesHouchesReader tmp_reader{filename};
-// reader().heprup.XSECUP = std::vector<double>{0};
- // while(tmp_reader.read_event()){
- // ++num_events_;
- // num_trials_ += std::stod(tmp_reader.hepeup().attributes.at("trials"));
- // reader().heprup.XSECUP.front() += tmp_reader.hepeup().XWGTUP;
- // }
- // reader().heprup.XSECUP.front() /= num_trials_;
Generator=get_generator(reader().heprup, reader().headerBlock);
int nprup=reader().heprup.NPRUP;
num_trials_.resize(nprup);
reader().heprup.XSECUP.resize(nprup);
// check if XSECUP is already set
if (reader().heprup.XSECUP.front()>0) {// if it is set
CalculateXSECUP=false;
}
- // // For IDWTUP == 1 or 4 we assume avg(weight)=xs
- // // With the modified weights we have in Sherpa sum(weight)=xs
- // // -> overwrite IDWTUP to "something neutral"
+ // For IDWTUP == 1 or 4 we assume avg(weight)=xs
+ // With the modified weights we have in Sherpa sum(weight)=xs
+ // -> overwrite IDWTUP to "something neutral"
reader().heprup.IDWTUP = reader().heprup.IDWTUP>0?3:-3;
}
bool HEJLHEReader::read_event() {
if(!LesHouchesReader::read_event()) return false;
-#if 0
- if (Generator==HEJ::generator::Sherpa) {
- ++num_events_;
- // if (HEJLHEReader::get_CalcXSECUP())
- reader().heprup.XSECUP.front() *= num_trials_;
- num_trials_ += std::stod(hepeup().attributes.at("trials"));
- // if (HEJLHEReader::get_CalcXSECUP()) {
- reader().heprup.XSECUP.front() += hepeup().XWGTUP;
- reader().heprup.XSECUP.front() /= num_trials_;
- // }
- } else {
- ++num_events_;
- if (HEJLHEReader::get_CalcXSECUP()) {
- reader().heprup.XSECUP.front() += hepeup().XWGTUP;
- }
- return true;
- }
-#else
int nprup=reader().heprup.NPRUP;
int position=std::max(1,std::min(hepeup().IDPRUP,nprup))-1; // 0..nprup-1
if (Generator == HEJ::generator::Sherpa) {
++num_events_;
- // if (HEJLHEReader::get_CalcXSECUP())
reader().heprup.XSECUP.at(position) *= num_trials_.at(position);
num_trials_.at(position) +=
std::stod(hepeup().attributes.at("trials"));
- // if (HEJLHEReader::get_CalcXSECUP()) {
reader().heprup.XSECUP.at(position) += hepeup().XWGTUP;
reader().heprup.XSECUP.at(position) /=
num_trials_.at(position);
- // }
} else {
++num_events_;
if (HEJLHEReader::get_CalcXSECUP()) {
reader().heprup.XSECUP.at(position) += hepeup().XWGTUP;
}
}
-#endif
- return true;
+ return true;
}
} // namespace HEJ
diff --git a/src/RivetAnalysis.cc b/src/RivetAnalysis.cc
index 40deecb..13e7ad3 100644
--- a/src/RivetAnalysis.cc
+++ b/src/RivetAnalysis.cc
@@ -1,218 +1,208 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "HEJ/RivetAnalysis.hh"
#include "HEJ/ConfigFlags.hh"
#include "HEJ/utility.hh"
#ifdef HEJ_BUILD_WITH_RIVET
#include <cstddef>
#include <ostream>
#include <utility>
#include "yaml-cpp/yaml.h"
#include "Rivet/AnalysisHandler.hh"
#include "Rivet/Config/RivetConfig.hh"
#ifdef RIVET_ENABLE_HEPMC_3
#include "HepMC3/GenEvent.h"
#include "HEJ/HepMC3Interface.hh"
#else // rivet with HepMC 2
#include "HEJ/HepMC2Interface.hh"
#include "HepMC/GenEvent.h"
#endif // RIVET_ENABLE_HEPMC_3
#include "HEJ/Event.hh"
#include "HEJ/exceptions.hh"
#endif
namespace HEJ {
std::unique_ptr<Analysis> RivetAnalysis::create(
YAML::Node const & config, LHEF::HEPRUP const & heprup
){
return std::make_unique<RivetAnalysis>(config, heprup);
}
}
#ifdef HEJ_BUILD_WITH_RIVET
namespace HEJ {
#ifdef RIVET_ENABLE_HEPMC_3
using HepMCInterface=HepMC3Interface;
#else
using HepMCInterface=HepMC2Interface;
#endif
struct RivetAnalysis::rivet_info {
rivet_info(std::unique_ptr<Rivet::AnalysisHandler> && h,
std::string && s, HepMCInterface && m):
handler{std::move(h)}, name{std::move(s)}, hepmc{std::move(m)}
{}
// AnalysisHandler doesn't allow a copy constructor -> use ptr
std::unique_ptr<Rivet::AnalysisHandler> handler;
std::string name;
HepMCInterface hepmc;
};
RivetAnalysis::RivetAnalysis(YAML::Node const & config,
LHEF::HEPRUP heprup
):
heprup_{std::move(heprup)},
first_event_(true)
{
// assert that only supported options are provided
if(!config.IsMap()) throw invalid_type{"Rivet config has to be a map!"};
for(auto const & entry: config){
const auto name = entry.first.as<std::string>();
if(name != "output" && name != "rivet")
throw unknown_option{"Unknown option \'"+name+"\' provided to rivet."};
}
// get output name
if(!config["output"])
throw std::invalid_argument{
"No output was provided to rivet. "
"Either give an output or deactivate rivet."
};
if(!config["output"].IsScalar())
throw invalid_type{
"Output for rivet must be a scalar."
};
output_name_ = config["output"].as<std::string>();
// read in analyses
auto const & name_node = config["rivet"];
switch(name_node.Type()){
case YAML::NodeType::Scalar:
analyses_names_.push_back(name_node.as<std::string>());
break;
case YAML::NodeType::Sequence:
for(YAML::const_iterator it = name_node.begin(); it != name_node.end(); ++it){
analyses_names_.push_back(it->as<std::string>());
}
break;
default:
throw std::invalid_argument{
"No analysis was provided to rivet. "
"Either give an analysis or deactivate rivet."
};
}
}
// it is a bit ugly that we can't do this directly in `initialise`
void RivetAnalysis::init(Event const & event){
#ifdef HEJ_USE_RIVET2
rivet_runs_.reserve(event.variations().size()+1);
#else
ignore(event); // shut up compiler
#endif
rivet_runs_.emplace_back( std::make_unique<Rivet::AnalysisHandler>(),
std::string(""), HepMCInterface(heprup_)
);
rivet_runs_.back().handler->addAnalyses(analyses_names_);
#ifdef HEJ_USE_RIVET2
//! scale variation in rivet 2 through multiple analyses
if( !event.variations().empty() ){
for(auto const & vari : event.variations()){
rivet_runs_.emplace_back( std::make_unique<Rivet::AnalysisHandler>(),
"."+to_simple_string(*vari.description), HepMCInterface(heprup_)
);
rivet_runs_.back().handler->addAnalyses(analyses_names_);
}
}
#endif
}
void RivetAnalysis::fill(Event const & event, Event const & /*unused*/){
if(first_event_){
first_event_=false;
init(event);
}
auto hepmc_event(rivet_runs_.front().hepmc(event));
rivet_runs_.front().handler->analyze(hepmc_event);
#ifdef HEJ_USE_RIVET2
for(std::size_t i = 1; i < rivet_runs_.size(); ++i){
auto & run = rivet_runs_[i];
run.hepmc.set_central(hepmc_event, event, i-1);
run.handler->analyze(hepmc_event);
}
#endif
}
void RivetAnalysis::scale(double xsscale){
- // used to scale the cross sections and histograms before .finalise in case the
- // normalisation is unknown until the end of the run
+ // used to scale the cross sections and histograms before .finalise in case the
+ // normalisation is unknown until the end of the run
for(auto & run: rivet_runs_){
- // get yodaOS:
+ // get yoda analysis objects:
std::vector<YODA::AnalysisObjectPtr> yodaAOS=run.handler->getYodaAOs(true);
for (auto& t : yodaAOS) {
- // print out name of yodaAO:
-// std::cout<<t->name()<<std::endl;//
-
auto * obj = dynamic_cast<YODA::Fillable*>(t.get());
if(obj) {
obj->scaleW(xsscale);
}
}
}
}
void RivetAnalysis::finalise(){
for(auto & run: rivet_runs_){
run.handler->finalize();
run.handler->writeData(output_name_+run.name+std::string(".yoda"));
}
}
} // namespace HEJ
#else // no rivet => create empty analysis
namespace Rivet {
class AnalysisHandler {};
}
namespace HEJ {
struct RivetAnalysis::rivet_info{};
RivetAnalysis::RivetAnalysis(
YAML::Node const & /*config*/, LHEF::HEPRUP /*heprup*/
){
ignore(first_event_);
throw std::invalid_argument(
"Failed to create RivetAnalysis: "
"HEJ 2 was built without rivet support"
);
}
void RivetAnalysis::init(Event const & /*event*/){}
void RivetAnalysis::fill(Event const & /*event*/, Event const & /*unused*/){}
- void RivetAnalysis::scale(double xsscale){
- // used to scale the cross sections and histograms before .finalise in case the
- // normalisation is unknown until the end of the run
- std::cout <<"RivetAnalysis Warning:: Function scale does nothing "<<xsscale<<std::endl;
- // for(auto & run: rivet_runs_){
- // Rivet::Scatter1DPtr temp=run.handler->crossSection();
- // }
- }
+ void RivetAnalysis::scale(double /*xsscale*/){}
void RivetAnalysis::finalise(){}
} // namespace HEJ
#endif
namespace HEJ {
RivetAnalysis::~RivetAnalysis() = default;
}