Page MenuHomeHEPForge
Files F19251732

test_unweighter.cc
No OneTemporary
Actions

Size
4 KB
Referenced Files
None
Subscribers
None

test_unweighter.cc
View Options

/**
* \authors The HEJ collaboration (see AUTHORS for details)
* \date 2019-2020
* \copyright GPLv2 or later
*/
#include "hej_test.hh"
#include <cmath>
#include <iostream>
#include <cstdlib>
#include <string>
#include <vector>
#include "HEJ/CrossSectionAccumulator.hh"
#include "HEJ/Event.hh"
#include "HEJ/EventReader.hh"
#include "HEJ/Mixmax.hh"
#include "HEJ/Unweighter.hh"
#include "fastjet/JetDefinition.hh"
namespace {
const fastjet::JetDefinition JET_DEF{fastjet::kt_algorithm, 0.4};
const double MIN_PT{30.};
const double SAMPLE_RATIO{10.};
const double MAX_DEV{2.};
bool compare_xs(
HEJ::XSWithError<double> const & xs1, HEJ::XSWithError<double> const & xs2
){
std::cout << xs1.value << "+/-" << xs1.error << " vs. "
<< xs2.value << "+/-" << xs2.error << std::endl;
return std::abs(xs1.value/xs2.value-1.)<xs1.error;
}
} // namespace
int main(int argc, char** argv) {
using std::size_t;
if(argc != 2) {
std::cerr << "Usage: " << argv[0] << " event_file\n";
return EXIT_FAILURE;
}
std::string file{argv[1]};
auto reader = HEJ::make_reader(file);
// number of events
auto nevents{reader->number_events()};
if(!nevents) {
auto t_reader = HEJ::make_reader(file);
nevents = 0;
while(t_reader->read_event()) ++(*nevents);
}
ASSERT(*nevents>SAMPLE_RATIO);
const size_t size_sample = std::floor(*nevents/SAMPLE_RATIO);
HEJ::Mixmax ran{};
// no unweighting
HEJ::CrossSectionAccumulator xs_base;
std::vector<HEJ::Event> all_evts;
// full unweighting
HEJ::CrossSectionAccumulator xs_max;
HEJ::Unweighter unw_max;
size_t n_max{0};
// midpoint on full sample
HEJ::CrossSectionAccumulator xs_mid;
HEJ::Unweighter unw_mid;
size_t n_mid{0};
// calc max from partial sample
HEJ::CrossSectionAccumulator xs_pmax;
HEJ::Unweighter unw_pmax;
size_t n_pmax{0};
// midpoint on partial sample
HEJ::CrossSectionAccumulator xs_pmid;
HEJ::Unweighter unw_pmid;
size_t n_pmid{0};
// initialise sample
for(size_t n = 0; n < size_sample; ++n){
if(!reader->read_event()){
std::cerr << "Sample size bigger than event sample\n";
return EXIT_FAILURE;
}
const HEJ::Event event{
HEJ::Event::EventData{reader->hepeup()}.cluster(JET_DEF, MIN_PT)
};
xs_base.fill(event);
all_evts.push_back(event);
}
// calculate partial settings
unw_pmax.set_cut_to_maxwt(all_evts);
unw_pmid.set_cut_to_peakwt(all_evts, MAX_DEV);
for(auto const & ev: unw_pmax.unweight(all_evts, ran)){
xs_pmax.fill(ev);
++n_pmax;
}
for(auto const & ev: unw_pmid.unweight(all_evts, ran)){
xs_pmid.fill(ev);
++n_pmid;
}
while(reader->read_event()){
const HEJ::Event event{
HEJ::Event::EventData{reader->hepeup()}.cluster(JET_DEF, MIN_PT)
};
xs_base.fill(event);
auto ev{ unw_pmid.unweight(event, ran) };
if(ev){
xs_pmid.fill(*ev);
++n_pmid;
}
ev = unw_pmax.unweight(event, ran);
if(ev){
xs_pmax.fill(*ev);
++n_pmax;
}
all_evts.push_back(event);
}
unw_max.set_cut_to_maxwt(all_evts);
unw_mid.set_cut_to_peakwt(all_evts, MAX_DEV);
for(auto const & ev: unw_max.unweight(all_evts, ran)){
// make sure all the events have the same weight
ASSERT( std::abs( std::abs(unw_max.get_cut()/ev.central().weight)-1. ) < 10e-16);
xs_max.fill(ev);
++n_max;
}
for(auto const & ev: unw_mid.unweight(all_evts, ran)){
xs_mid.fill(ev);
++n_mid;
}
// sanity check number of events
ASSERT( !all_evts.empty());
ASSERT( n_pmax > 0);
ASSERT( n_max > 0);
ASSERT( n_pmid > 0);
ASSERT( n_mid > 0);
ASSERT( n_pmax < all_evts.size() );
ASSERT( n_max < n_pmax );
ASSERT( n_pmid < all_evts.size() );
ASSERT( n_mid < all_evts.size() );
ASSERT( n_max < n_mid );
std::cout << "all_evts.size() " << all_evts.size() << " n_pmax " << n_pmax <<
" n_max " << n_max << " n_pmid " << n_pmid << " n_mid " << n_mid << std::endl;
// control xs (in circle)
ASSERT(compare_xs( xs_base.total(), xs_pmax.total() ));
ASSERT(compare_xs( xs_pmax.total(), xs_max.total() ));
ASSERT(compare_xs( xs_max.total() , xs_pmid.total() ));
ASSERT(compare_xs( xs_pmid.total(), xs_mid.total() ));
ASSERT(compare_xs( xs_mid.total() , xs_base.total() ));
return EXIT_SUCCESS;
}

File Metadata

Mime Type
text/x-c
Expires
Tue, Sep 30, 6:09 AM (1 d, 6 h)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
6525286
Default Alt Text
test_unweighter.cc (4 KB)

Event Timeline

HEPForge.org · .