Page MenuHomeHEPForge
Files F9501778

No OneTemporary
Actions

View Options

diff --git a/src/LesHouchesWriter.cc b/src/LesHouchesWriter.cc
index 2fae860..0b2693e 100644
--- a/src/LesHouchesWriter.cc
+++ b/src/LesHouchesWriter.cc
@@ -1,123 +1,124 @@
/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019
* \copyright GPLv2 or later
*/
#include <cassert>
#include <utility>
#include <vector>
#include "HEJ/Event.hh"
#include "HEJ/event_types.hh"
#include "HEJ/LesHouchesWriter.hh"
#include "HEJ/utility.hh"
namespace HEJ{
namespace{
template<class T, class... Args>
std::unique_ptr<T> make_unique(Args&&... a){
return std::unique_ptr<T>{new T{std::forward<Args>(a)...}};
}
size_t to_index(event_type::EventType const type){
return type==0?0:floor(log2(type))+1;
}
}
LesHouchesWriter::LesHouchesWriter(
std::string const & file, LHEF::HEPRUP heprup
):
out_{file, std::fstream::in | std::fstream::out | std::fstream::trunc},
writer_{HEJ::make_unique<LHEF::Writer>(out_)}
{
if(! out_.is_open()){
throw std::ios_base::failure("Failed to open " + file);
};
// scientific style is needed to allow rewriting the init block
out_ << std::scientific;
writer_->heprup = std::move(heprup);
// lhe Standard: IDWTUP (negative => weights = +/-)
// IDWTUP: HEJ -> SHG/Pythia/next program
// 1: weighted->unweighted, xs = mean(weight), XMAXUP given
// 2: weighted->unweighted, xs = XSECUP, XMAXUP given
// 3: unweighted (weight=+1)->unweighted, no additional information
// 4: weighted->weighted, xs = mean(weight)
//
// None of these codes actually match what we want:
// 1 and 4 require xs = mean(weight), which is impossible until after generation
// 2 tells the SHG to unweight our events, which is wasteful
// 3 claims we produce unweighted events, which is both wasteful _and_
// impossible until after generation (we don't know the maximum weight before)
//
- // For the time being, we choose 3. If the consumer (like Pythia) assumes
+ // For the time being, we choose -3. If the consumer (like Pythia) assumes
// weight=+1, the final weights have to be corrected by multiplying with
- // the original weight we provided.
+ // the original weight we provided. We are also often use NLO-PDFs which can
+ // give negative weights, hence the native IDWTUP.
//
- writer_->heprup.IDWTUP = 3;
+ writer_->heprup.IDWTUP = -3;
const int max_number_types = to_index(event_type::last_type)+1;
writer_->heprup.NPRUP = max_number_types;
// ids of event types
writer_->heprup.LPRUP.clear();
writer_->heprup.LPRUP.reserve(max_number_types);
writer_->heprup.LPRUP.emplace_back(0);
for(size_t i=event_type::first_type+1; i<=event_type::last_type; i*=2)
writer_->heprup.LPRUP.emplace_back(i);
// use placeholders for unknown init block values
// we can overwrite them after processing all events
writer_->heprup.XSECUP = std::vector<double>(max_number_types, 0.);
writer_->heprup.XERRUP = std::vector<double>(max_number_types, 0.);
writer_->heprup.XMAXUP = std::vector<double>(max_number_types, 0.);
write_init();
}
void LesHouchesWriter::write(Event const & ev){
assert(writer_ && out_.is_open());
const double wt = ev.central().weight;
writer_->hepeup = HEJ::to_HEPEUP(std::move(ev), &heprup());
writer_->writeEvent();
assert(heprup().XSECUP.size() > to_index(ev.type()));
heprup().XSECUP[to_index(ev.type())] += wt;
heprup().XERRUP[to_index(ev.type())] += wt*wt;
if(wt > heprup().XMAXUP[to_index(ev.type())]){
heprup().XMAXUP[to_index(ev.type())] = wt;
}
}
// this function is called after overwritting the Les Houches init block
// assert that we have overwritten *exactly* the init block,
// i.e. we are at the end of the file or an intact event block is next
void assert_next_event_intact(std::istream & out){
(void) out; // suppress compiler warnings if not in debug mode
#ifndef NDEBUG
std::string line;
getline(out, line);
assert(out.eof() || line.rfind("<event", 0) == 0);
#endif
}
void LesHouchesWriter::rewrite_init(){
assert(writer_ && out_.is_open());
// replace placeholder entries
const auto pos = out_.tellp();
out_.seekp(0);
write_init();
assert_next_event_intact(out_);
out_.seekp(pos);
}
LesHouchesWriter::~LesHouchesWriter(){
assert(writer_ && out_.is_open());
for(auto & xs_err: heprup().XERRUP)
{
xs_err = sqrt(xs_err);
}
rewrite_init();
}
}

File Metadata

Mime Type
text/x-diff
Expires
Sun, Feb 23, 3:03 PM (19 m, 24 s)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4486822
Default Alt Text
(4 KB)

Event Timeline