Page MenuHomeHEPForge
Files F8725694

No OneTemporary
Actions

View Options

diff --git a/analysis-plugins/include/StandardAnalysis_impl.hh b/analysis-plugins/include/StandardAnalysis_impl.hh
index 3a0b19b..7017619 100644
--- a/analysis-plugins/include/StandardAnalysis_impl.hh
+++ b/analysis-plugins/include/StandardAnalysis_impl.hh
@@ -1,525 +1,525 @@
#pragma once
#include <iostream>
#include <numeric>
#include "StandardAnalysis.hh"
#include "ROOTWriter.hh"
#include "ROOTReader.hh"
#include "TDirectory.h"
#include "TKey.h"
#include "TCanvas.h"
#include "TImage.h"
#include "TFile.h"
#include "yaml-cpp/yaml.h"
namespace RHEJ_analyses{
namespace detail{
inline
bool disable_TH1_reference(){
TH1::AddDirectory(kFALSE);
std::cout << "Disabled central ROOT histogram references\n";
return true;
};
}
using histogram = StandardAnalysis::histogram_type;
inline
StandardAnalysis::hist_pack::hist_pack(
const char * title, const std::string & id,
int numbins, double min, double max
):
h((title + id).c_str(), title, numbins, min, max),
hwt(
(title + ("wt2" + id)).c_str(),
(title + std::string("wt2")).c_str(),
numbins, min, max
),
hN(
(title + ("N" + id)).c_str(),
(title + std::string("N")).c_str(),
numbins, min, max
)
{
static const bool TH1_reference_disabled = detail::disable_TH1_reference();
(void) TH1_reference_disabled; // avoid compiler warnings
h.sumw2();
}
inline
void StandardAnalysis::hist_pack::write() const{
h.write();
hwt.write();
hN.write();
}
inline
int StandardAnalysis::hist_pack::nbins() const{
return h.nbins();
}
inline
void StandardAnalysis::hist_pack::fill(double val, double wt){
const int binnumber = h.x_axis()->FindBin(val);
const double binwidth = h.x_axis()->GetBinWidth(binnumber);
h.fill(val, wt/binwidth);
hwt.fill(val, wt*wt/binwidth);
hN.fill(val, 1);
}
inline
StandardAnalysis::EventHists::EventHists():
EventHists{""} {}
StandardAnalysis::EventHists::EventHists(std::string const & id){
#define NEWHIST(X,UUID,BINS,MIN,MAX) hist_[X] = hist_pack{#X, UUID, BINS, MIN, MAX};
NEWHIST(total,id,1,-.5,.5);
NEWHIST(np,id,35,-.5,34.5);
NEWHIST(njets,id,7,.5,7.5);
NEWHIST(pt1,id,100,0.,1000.);
NEWHIST(pt2,id,100,0.,1000.);
NEWHIST(pth1,id,100,0.,1000.);
NEWHIST(pth2,id,100,0.,1000.);
NEWHIST(qt,id,100,0.,99.);
NEWHIST(ya,id,40,-5.,5.);
NEWHIST(yb,id,40,-5.,5.);
NEWHIST(ydif,id,32,0.,8.);
NEWHIST(yEW,id,40,-5.,5.);
NEWHIST(ptEW,id,20,0.,200.);
#undef NEWHIST
hist_wt_ = histogram{("wt" + id).c_str(),"wt",800,-30.,50.};
hist_wtwt_ = histogram{("wtwt" + id).c_str(),"wtwt",800,-30.,50.};
}
inline
void StandardAnalysis::EventHists::fill(KinInfo const & i, double weight){
const double ymin = i.jets_y.front().rapidity();
const double ymax = i.jets_y.back().rapidity();
assert(ymax > ymin);
const auto EW_boson = std::find_if(
begin(i.out), end(i.out),
[](RHEJ::Sparticle const & s){ return is_AWZH_boson(s); }
);
const size_t npartons = (EW_boson==end(i.out))?
i.out.size():
i.out.size()-1
;
hist_[total].fill(0., weight);
hist_[np].fill(npartons, weight);
for(size_t jets = 2; jets <= i.jets_y.size(); ++jets){
hist_[njets].fill(jets, weight);
}
hist_[pt1].fill(i.jets_y.front().perp(), weight);
hist_[pt2].fill(i.jets_y.back().perp(), weight);
hist_[pth1].fill(i.jets_pt.front().perp(), weight);
hist_[pth2].fill(i.jets_pt[1].perp(), weight);
hist_[qt].fill(i.q.perp(), weight);
hist_[ya].fill(ymin, weight);
hist_[yb].fill(ymax, weight);
hist_[ydif].fill(ymax - ymin, weight);
if(EW_boson!=end(i.out)){
hist_[yEW].fill(EW_boson->rapidity(), weight);
hist_[ptEW].fill(EW_boson->perp(), weight);
}
const double awt = std::abs(weight);
hist_wt_.fill(std::log(awt), 1);
hist_wtwt_.fill(std::log(awt), awt);
}
inline
void StandardAnalysis::EventHists::write() const{
for(auto hist: hist_){
hist.write();
}
hist_wt_.write();
hist_wtwt_.write();
}
namespace detail{
inline
double get_log_wt_cut(YAML::Node const & cut){
return cut?
cut.as<double>():
std::numeric_limits<double>::infinity();
}
}
inline
StandardAnalysis::StandardAnalysis(YAML::Node const & parameters):
log_wt_cut_{detail::get_log_wt_cut(parameters["log wt cut"])},
output_{parameters["output"].as<std::string>()}
{}
inline
StandardAnalysis::StandardAnalysis(
std::string const & filename
){
restore(filename);
}
inline
StandardAnalysis::~StandardAnalysis(){
if(!output_.empty()) store();
}
inline
StandardAnalysis::StandardAnalysis(StandardAnalysis const & other){
*this = other;
}
inline
StandardAnalysis::StandardAnalysis(StandardAnalysis && other){
*this = std::move(other);
}
inline
StandardAnalysis & StandardAnalysis::operator=(StandardAnalysis const & other){
if(this == &other) return *this;
central_ = other.central_;
hist_ = other.hist_;
log_wt_cut_ = other.log_wt_cut_;
output_ = other.output_;
bad_events_ = other.bad_events_;
orig_central_ = other.orig_central_;
orig_hist_ = other.orig_hist_;
orig_log_wt_cut_ = other.orig_log_wt_cut_;
update_hist_by_name();
return *this;
}
StandardAnalysis & StandardAnalysis::operator=(StandardAnalysis && other){
if(this == &other) return *this;
central_ = std::move(other.central_);
hist_ = std::move(other.hist_);
log_wt_cut_ = std::move(other.log_wt_cut_);
output_ = std::move(other.output_);
bad_events_ = std::move(other.bad_events_);
orig_central_ = std::move(other.orig_central_);
orig_hist_ = std::move(other.orig_hist_);
orig_log_wt_cut_ = std::move(other.orig_log_wt_cut_);
update_hist_by_name();
return *this;
}
inline
std::vector<std::string> StandardAnalysis::histogram_names() const{
std::vector<std::string> result;
for(auto const & entry: hist_by_name_) result.emplace_back(entry.first);
return result;
};
inline
std::vector<RHEJ::DynamicAnalysis::Parameter> StandardAnalysis::parameters() const{
double min = log_wt_cut_;
double max = log_wt_cut_;
if(! bad_events_.empty()){
min = std::min(min, std::log(std::abs(begin(bad_events_)->central().weight)));
- max = std::max(max, std::log(std::abs(rbegin(bad_events_)->central().weight)));
+ max = std::max(max, std::log(std::abs(bad_events_.rbegin()->central().weight)));
}
return {RHEJ::DynamicAnalysis::Parameter{"log wt cut", log_wt_cut_, min, max}};
}
inline
void StandardAnalysis::draw_histogram(
std::string const & name, std::string const & filename
) const{
const auto hist = hist_by_name_.find(name);
if(hist == end(hist_by_name_)){
throw std::invalid_argument{
"histogram with title "+ name + " does not exist"
};
}
TCanvas c;
hist->second.get().draw();
const auto img = std::unique_ptr<TImage>{TImage::Create()};
img->FromPad(&c);
img->WriteImage(filename.c_str(), TImage::EImageFileTypes::kPng);
}
inline
bool StandardAnalysis::update_parameter(std::string const & name, double value){
if(name != "log wt cut"){
throw std::invalid_argument{"invalid parameter " + name};
}
Event search_event;
search_event.central().weight = std::exp(log_wt_cut_);
const auto old_end = bad_events_.lower_bound(search_event);
search_event.central().weight = std::exp(value);
const auto new_end = bad_events_.lower_bound(search_event);
if(new_end == old_end){
log_wt_cut_ = value;
return false; // no update required
}
if(value < log_wt_cut_){
// we would have to remove events from histograms
// instead, reset and add all events again
reset();
update_parameter(name, value);
return true;
}
log_wt_cut_ = value;
for(auto it = old_end; it != new_end; ++it) fill(*it);
return true;
}
inline
void StandardAnalysis::fill(Event const & ev){
if(std::log(std::abs(ev.central().weight)) > log_wt_cut_){
bad_events_.emplace(ev);
return;
}
const std::vector<fastjet::PseudoJet> partons =
to_PseudoJet(filter_partons(ev.outgoing()));
auto jets_y = sorted_by_rapidity(ev.jets());
auto q =
std::accumulate(partons.begin(), partons.end(), fastjet::PseudoJet{})
- std::accumulate(jets_y.begin(), jets_y.end(), fastjet::PseudoJet{})
;
const KinInfo info = {
ev.outgoing(),
std::move(jets_y),
sorted_by_pt(ev.jets()),
std::move(q)
};
assert(info.jets_y.size() == info.jets_pt.size());
assert(info.jets_y.size() >= 2);
central_.fill(info, ev.central().weight);
if(hist_.empty()){
hist_.reserve(ev.variations().size());
for(size_t i = 0; i < ev.variations().size(); ++i){
hist_.emplace_back('_' + std::to_string(i));
}
}
if(hist_.size() != ev.variations().size()){
throw std::invalid_argument(
"number of weights " + std::to_string(ev.variations().size())
+ " does not match number of histograms "
+ std::to_string(hist_.size())
);
}
for(size_t i = 0; i < ev.variations().size(); ++i){
hist_[i].fill(info, ev.variations(i).weight);
}
}
inline
void StandardAnalysis::write_settings(TFile &){
TTree settings{"settings", "settings"};
settings.Branch("log_wt_cut", &log_wt_cut_, "wtwt_cut/D");
settings.Fill();
settings.Write();
}
inline
void StandardAnalysis::read_settings(TFile & in){
auto tree = static_cast<TTree*>(in.Get("settings"));
if(tree == nullptr){
throw std::runtime_error("no TTree with name \"settings\"");
}
tree->SetBranchAddress("log_wt_cut", &log_wt_cut_);
if(tree->GetEntries() != 1){
throw std::runtime_error("wrong number of entries in \"settings\"");
}
tree->GetEntry(0);
}
inline
void StandardAnalysis::write_histograms(TFile & out){
const auto dir = out.mkdir("histograms");
assert(dir != nullptr);
dir->cd();
central_.write();
for(auto const & hist: hist_) hist.write();
out.cd();
}
inline
void StandardAnalysis::write_events(TFile &){
ROOTWriter writer{"events"};
for(auto && event: bad_events_) writer.write(event);
}
namespace detail{
struct weight_less{
bool operator()(Event const & a, Event const & b){
return a.central().weight < b.central().weight;
}
bool operator()(Event const & ev, double wt){
return ev.central().weight < wt;
}
bool operator()(double wt, RHEJ::Event const & ev){
return wt < ev.central().weight;
}
};
}
inline
void StandardAnalysis::read_events(TFile & in){
assert(in.IsOpen());
auto event_tree = static_cast<TTree*>(in.Get("events"));
if(event_tree == nullptr){
throw std::runtime_error("no TTree with name \"events\"");
}
ROOTReader event_reader{event_tree};
bad_events_ = std::set<Event, detail::AbsWtLess>(
std::begin(event_reader), std::end(event_reader)
);
}
inline
void StandardAnalysis::store(){
TFile out{output_.c_str(), "RECREATE"};
write_settings(out);
write_histograms(out);
write_events(out);
}
inline
void StandardAnalysis::restore(std::string const & filename){
TFile in{filename.c_str()};
read_settings(in);
read_histograms(in);
read_events(in);
orig_central_ = central_;
orig_hist_ = hist_;
orig_log_wt_cut_ = log_wt_cut_;
}
namespace detail{
using histogram_type = StandardAnalysis::histogram_type;
inline
size_t num_hist_indices(
std::map<std::string, Histogram> const & hists,
std::string const & name
){
size_t result = 0;
while(hists.find(name + "_" + std::to_string(result)) != hists.end()){
++result;
}
return result;
}
inline
std::map<std::string, histogram_type> read_hist_map(TDirectory & in){
TIter keys = in.GetListOfKeys();
TKey* key;
std::map<std::string, histogram_type> hists;
while( (key = static_cast<TKey*>(keys())) ){
auto new_hist = dynamic_cast<TH1D*>(key->ReadObj());
if(! new_hist){
throw std::runtime_error("Encountered unexpected object in ROOT file");
}
hists.emplace(new_hist->GetName(), histogram_type{new_hist});
}
return hists;
}
}
inline
void StandardAnalysis::read_histograms(TFile & in){
const auto dir = in.GetDirectory("histograms");
if(dir == nullptr){
throw std::runtime_error{"no directory \"histograms\""};
}
restore_histograms(detail::read_hist_map(*dir));
update_hist_by_name();
}
inline
void StandardAnalysis::restore_histograms(
std::map<std::string, histogram_type> && hists
){
hist_.resize(detail::num_hist_indices(hists, hist_names[first_histogram]));
for(size_t id = first_histogram; id <= last_histpack; ++id){
central_.hist_[id].h = std::move(hists[hist_names[id]]);
for(size_t i = 0; i < hist_.size(); ++i){
const std::string name = hist_names[id] + ("_" + std::to_string(i));
hist_[i].hist_[id].h = std::move(hists[name]);
}
std::string name = hist_names[id] + std::string{"wt2"};
central_.hist_[id].hwt = std::move(hists[name]);
for(size_t i = 0; i < hist_.size(); ++i){
const std::string name = hist_names[id] + ("wt2_" + std::to_string(i));
hist_[i].hist_[id].hwt = std::move(hists[name]);
}
name = hist_names[id] + std::string{"N"};
central_.hist_[id].hN = std::move(hists[name]);
for(size_t i = 0; i < hist_.size(); ++i){
const std::string name = hist_names[id] + ("N_" + std::to_string(i));
hist_[i].hist_[id].hN = std::move(hists[name]);
}
}
central_.hist_wt_ = std::move(hists["wt"]);
central_.hist_wtwt_ = std::move(hists["wtwt"]);
for(size_t i = 0; i < hist_.size(); ++i){
std::string name = "wt_" + std::to_string(i);
hist_[i].hist_wt_ = std::move(hists[name]);
name = "wtwt_" + std::to_string(i);
hist_[i].hist_wtwt_ = std::move(hists[name]);
}
}
inline
void StandardAnalysis::update_hist_by_name(){
hist_by_name_.clear();
add_to_hist_by_name(central_);
for(auto const & hist: hist_) add_to_hist_by_name(hist);
}
inline
void StandardAnalysis::add_to_hist_by_name(EventHists const & hist){
for(auto const & h: hist.hist_) add_to_hist_by_name(h);
hist_by_name_.emplace(hist.hist_wt_.title(), Hist_ref{hist.hist_wt_});
hist_by_name_.emplace(hist.hist_wtwt_.title(), Hist_ref{hist.hist_wtwt_});
}
inline
void StandardAnalysis::add_to_hist_by_name(hist_pack const & pack){
hist_by_name_.emplace(pack.h.title(), Hist_ref{pack.h});
hist_by_name_.emplace(pack.hwt.title(), Hist_ref{pack.hwt});
hist_by_name_.emplace(pack.hN.title(), Hist_ref{pack.hN});
}
inline
std::unique_ptr<RHEJ::Analysis> StandardAnalysis::create(
YAML::Node const & parameters
){
return std::unique_ptr<Analysis>{new StandardAnalysis{parameters}};
}
inline
std::unique_ptr<RHEJ::DynamicAnalysis> StandardAnalysis::create(
std::string const & filename
){
return std::unique_ptr<RHEJ::DynamicAnalysis>{
new StandardAnalysis{filename}
};
}
inline
void StandardAnalysis::reset(){
central_ = orig_central_;
hist_ = orig_hist_;
log_wt_cut_ = orig_log_wt_cut_;
update_hist_by_name();
}
constexpr
decltype(StandardAnalysis::hist_names) StandardAnalysis::hist_names;
}

File Metadata

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

Event Timeline